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Abstract 

In  this  investigation,  a  generalized  computer  code  for  the  simula¬ 
tion  of  fluid  transients  was  assembled,  verified,  and  applied.  The 
solution  method  developed  for  the  program  was  based  on  a  method  of 
characteristics  solution  of  the  equations  of  motion  for  one-dimensional 
fluid  flow  in  pipe  line  Sj^The  differential  equations  were  solved  using  a 
first-order  finite  difference  technique.  Boundary  conditions  for  the 
equations  of  motion  were  developed  or  directly  included  from  available 
component  models.  The  computer  routines  were  verified  by  comparison  of 
results  with  published  results  from  several  sources.  The  agreement  be¬ 
tween  the  results  of  this  study  and  the  published  data  was  quite  good  for 
a  wide  range  of  boundary  conditions  and  pipe  systems.  Recommendations 
were  made  for  improvement  of  some  models  and  the  replacement  of  others 
to  further  improve  the  accuracy  of  the  results. 
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DYNAMICS  OF  PROPELLANT  FEEDLINE  SYSTEMS 


I.  Introduction 

Fluid  systems  in  use  by  the  aerospace  industry  are  often  charac¬ 
terized  by  high  volume  flow  rates  or  high  system  pressures.  Furthermore, 
these  operating  conditions  are  coupled  with  saall  permissible  ranges  of 
variation  from  the  design  flow  and  pressure  operating  points  (1:1).  The 
net  result  of  the  two  factors  is  a  large  potential  for  performance  degra¬ 
dation  due  to  fluid  tranaients.  Within  the  space  prograa,  a  sometimes 
spectacular  manifestation  of  thia  phenomena  is  combustion  instability,  a 
performance  problem  frequently  caused  by  wave  motion  occuring  within  the 
propellant  feed  system  of  liquid  rocket  engines  (2:1). 

Previous  Work.  Many  attempts  have  been  made  to  model  the  dynamics 
of  liquid  propellant  feedlines.  All  of  the  models  developed  to  date, 
however,  can  be  catagorized  either  as  a  linearized  network  model  or  as  a 
time-step  simulation.  The  network  models  achieve  a  solution  by 
linearizing  the  equations  of  motion,  then  treating  the  transient  as  a 
perturbation  on  some  mean  flow.  Time-step  simulations  on  the  other  hand 
employ  nonlinear  finite  difference  forma  of  the  equations  of  motion, 
usually  obtained  by  the  method  of  characteristics  (3:1). 

One  of  the  first  to  apply  linearized  flow  equations  to  feedlines  was 
Rubin  (4).  He  employed  both  fluid  dynamic  theories  and  electrical  cir¬ 
cuit  analogies  to  achieve  a  solution.  Ryan  (3)  made  use  of  a  linearized 
model  to  investigate  structural  instabilities  coupled  to  fluid  transients 
in  the  S-IIc  engines  of  the  Saturn  V.  About  that  same  time  Johnson  (6) 
attempted  to  produce  a  generalized  computer  program  based  on  linear 


equations  to  predict  longitudinal  instability  in  propulsion  systems. 

More  recently,  Holster  and  Astleford  (7)  developed  a  general  analytical 
model  based  on  the  liner ized  flow  equations  developed  by  O' Souza  and 

i 

Oldenburger  (8).  i 

As  for  time-step  simulations.  Woods  (9)  applied  the  method  of 
characteristics  to  a  frictionless  model  to  study  fluid  fluctuations  in 
feedlines.  Wood  et  al .  (1)  developed  a  more  general  nonlinear  distri¬ 
buted  parameter  model  which  the  authors  named  the  "wave-plan"  method.  ! 

The  "wave-plan"  method  was  developed  to  predict  unsteady  flow  in  liquid 
filled  lines  and  was  later  used  to  create  a  generalized  digital  computer 
program  (10)  to  analyze  fluid  transients  in  liquid  rocket  feedlines. 

I 

Fashbaugh  and  Streeter  (11)  made  use  of  the  method  of  characteristics  to 
develop  their  own  digital  computer  program  to  investigate  transients  in 
propellant  feedlines  of  the  Titan  II  missile.  However,  one  of  the  most 
recent  digital  computer  programs  that  made  use  of  the  method  of  charac¬ 
teristics  did  not  arise  out  of  a  study  of  propellant  feedlines.  During 
the  Aircraft  Hydraulic  Systems  Dynamic  Analysis  Project,  the  McDonnell 
Aircraft  Company  developed  four  digital  computer  programs  for  simulating 
a  number  of  different  aspects  of  the  dynamics  of  aircraft  hydraulic 
systems.  One  of  the  four  programs,  named  Hydraulic  Transients  (HYTRAN) 

(12),  simulates  and  predicts  the  dynamic  response  of  a  hydraulic  system 
to  sudden  changes  in  load  flow  demand.  HYTRAN,  however,  is  a  very 
complex  computer  program,  where  the  user  must  be  very  familiar  with  tran¬ 
sient  phenomena  in  hydraulic  systems  and  know  how  to  interpret  the 
results. 

Objectives^  The  time-step  simulations  which  have  been  encountered 
in  the  literature,  so  far,  tended  to  be  either  very  complex  or  overly 


2 


restrictive.  Therefore,  the  objective  of  this  investigation  was  to 
assemble,  verify,  and  apply  a  simple  but  flexible  computer  program  for 
the  analysis  of  fluid  transients  in  pipelines.  Routines  in  the  program 
were  adapted  from  existing  code  whenever  possible.  In  addition,  this 
investigation  was  to  provide  a  source  from  which  future  investigators 
could  draw  references  on  work  of  a  similar  nature. 

Approach .  The  approach  used  in  this  investigation  was  not  comple¬ 
tely  new,  rather,  it  was  a  variation  on  the  method  used  in  several  of 
the  works  presented  earliear.  The  investigation  conducted  for  this  the¬ 
sis  proceeded  through  three  main  steps.  Step  one  was  the  derivation  of 
the  equations  that  would  form  the  basis  of  the  computer  program.  First 
the  finite  difference  forma  of  the  equations  of  motion  were  derived  by 
the  method  of  characteristics  (13,  14).  Then,  a  number  of  boundary  con¬ 
ditions  which  commonly  occur  in  propellant  feedlines  were  developed  to 
provide  for  a  complete  solution.  The  second  step  of  this  study  was  to 
present  the  systems  to  be  modeled.  It  was  at  this  point,  that  the  input 
data  for  each  system  was  described  and  the  appropriate  boundary  con¬ 
ditions  assigned.  And  finally,  the  last  step  was  to  analyze  the  results 
obtained  by  running  computer  simulations  which  made  use  of  the  systems 
modeled  in  step  two.  The  results  of  the  simulations  were  compared  with 
published  data  when  possible. 


II.  Method  of  Characteristics 


i 

In  this  section  a  numerical  solution  of  the  equations  that  were  ! 

used  to  model  unsteady  flow  in  pipelines  is  developed  by  the  method  of  l 

characteristics.  This  technique  was  used  to  transform  partial  differen¬ 
tial  equations  which  had  no  general  solution  into  particular  total  dif¬ 
ferential  equations  that  did.  The  equations  that  resulted  were  then  1 

integrated  to  provide  finite  difference  equations  by  which  the  pressure 
and  flow  velocity  within  a  pipeline  could  be  obtained  numerically.  The  I 

I 

process  by  which  these  equations  were  derived  was  not  new,  however  the 

final  form  of  the  finite  difference  equations  used  for  the  calculations  ! 

I 

was  different.  Unlike  previous  derivations,  the  final  form  was  in  terms  I 

of  pressure  and  volume  flow  rate,  not  head  or  pressure  and  velocity. 

Equations  of  Motion 

In  this  study,  transient  fluid  flow  was  represented  by  a  one¬ 
dimensional  model  with  time,  t,  and  distance  along  the  pipe  axis,  x,  as 
independent  variables,  and  pressure,  P,  and  mean  sectional  velocity,  V, 
as  dependent  varibles.  The  partial  differential  equations  may  be  derived 
by  the  application  of  the  momentum  and  of  the  continuity  principles  to  a 
constant  area  section  of  pipe  having  a  length,  dx  (13:694).  The  deriva¬ 
tion  which  follows  is  similar  to  derivations  presented  by  Swaffield  (13), 

Wylie  and  Streeter  (14),  and  Watters  (15). 

The  eqations  of  motion:- 

Li-Px/p  +gs ina+Wx+Vt+fV  |  V |  /2D-0  (2.1) 

and  continuity: 

L2-Pt+pa2Vx+VPx-0  (2.2) 

are  the  equations  which  will  be  used  in  this  development. 


I 


The  method  of  characteristics  proceeds  by  making  a  linear  com¬ 


bination  of  Eq.  (2.1)  and  (2.2)  using  an  unknown  multiplier,  X. 


L-L1  +  XL2-[Px/p+gsina+Wx+fV  V  /2D] +x[Pt+pa2Vx+VPx]-0 


(2.3) 


Regrouping  terms , 


X [ Px ( V+ 1 /x P ) +Pt 1 + [ Vx ( V + X Pa2 ) +V t ] +gs ina+ fV | V | / 2D-0 


(2.4) 


Eq.  (2.4)  can  be  simplified  by  the  appropriate  selection  of  the  two 


particular  values  of  X.  Since  P  and  V  are  functions  of  x  and  t,  from 


calculus 


dP/dt*Pxdx/dt+Pt 


dV/dt-Vxdx/dt+Vt 


(2.5) 


Upon  examination  of  equations  (2.4)  and  (2.5), 


dx/dt*V+l/XP“V+Xpa2 


(2.6) 


Thus,  Eq.  (2.4)  becomes  the  ordinary  differential  equation 


X dP/dt +dV/ dt +gs ina+fV | V | / 2D-0 


(2.7) 


Eq.  (2.6)  can  be  solved  to  obtain  the  two  particular  values  of  X, 


X  ■♦1/pi 


(2.8) 


Substituting  Eq.  (2.8)  back  into  Eq.  (2.6)  produces 


dx/dt-Vt-a 


(2.9) 


To  complete  the  transformation,  Eq.  (2.8)  is  substituted  into  Eq. 


(2.7).  The  characteristic  equations  which  result  are: 


dV/dt+(l/pa)dP/dt+gsina+fv|v|/2D-0 


(2.10) 


dx/ dt*V+a 


(2.11) 


dV/dt-(l/pa)dP/dt+gsina+fv|v|/2D-0 


(2.12) 


dx/dt“V-a 


(2.13) 


where  Eq.  (2.10)  and  (2.12)  are  valid  only  when  the  respective  equations, 


Eq.  (2.11)  and  (2.13)  are  valid. 


Finite  Difference  Equations 


Due  to  the  nonlinear  nature  of  the  characteristic  equations,  a 
finite  difference  approach  is  used  to  obtain  a  numerical  solution.  To 
this  point  the  derivation  has  paralleled  the  development  presented  by 
Watters  (15).  However,  Watters  completed  his  derivation  after  replacing 
pressure  by  piezometric  head,  P*pg(H-z).  Wylie  and  Streeter  (14)  make 
the  same  substitution  much  earlier  in  their  derivation.  This  develop¬ 
ment,  though,  will  retain  the  pressure  throughout,  in  a  manner  similar  to 
the  development  by  Swaf field  (13). 

By  multiplying  Eq.  (2.10)  through  (2.13)  by  dt ,  and  integrating  the 
equations  using  a  first-order  approximation,  the  finite  difference  forms 
are: 


Vp-VR+(Pp-PR)/pa+gsina(  tp-tR)+fVR|  VR|  ( tp-tg)/2D»0 

(2.14) 

xp-xR-(VR+a)(tp-tR) 

(2.15) 

Vp-Vs~(Pp-Ps)/ba+gsina( tp-tg)+fVg | Vg J ( tp-tg)/2D*0 

(2.16) 

Xp-Xg*(Vg-a) ( tp-tg) 

(2.17) 

Solving  for  Vp  and  writing  Eq.  (2.14)  and  (2.16)  in  terms 

of  volume 

flow  rate  in  place  of  velocity 

Qp*Cp-BPp 

(2.18) 

QpMCM+BPp 

(2.19) 

In  which 

cp«QR-BPR-RR-FFQR|QRj 

(2.20) 

cm"Qs_bps"rs”ffQs  1  Qs  1 

(2.21) 

where  B“A/pa,  and  RagAAtsina,  and  FF"fAt/2DA. 

Fig.  2-1,  which  graphically  represents  Eq.  (2.15)  and  (2.17),  shows 
the  characteristics  in  the  x-t  plane  along  which  Eq.  (2.14)  and  (2.16) 
are  integrated. 
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Fig.  2-1.  x-t  grid  £or  a  pipe  divided  into  N  reaches. 


To  determine  Che  mesh  of  the  grid  Che  method  of  specified  cime 
intervals  is  used.  In  order  Co  keep  Che  soluCion  numerically  stable,  Che 


CouranC  condicion  (14:58)  aueC  be  sacisfied 

&C(V+a)<Ax 


(2.22) 


The  aeChod  of  specified  cime  intervals  assumes  the  condicions  are  known 
aC  A,  B,  and  C  (Fig.  2-1)  eicher  from  Che  previous  Cime  seep  or  from  Che 
sCeady  solution.  Following  Che  usual  approach,  a  linear  inCerpolaCion  is 
Chen  made  Co  find  Che  pressure  and  Che  volume  flow  race  aC  points  R  and 
S.  The  stability  criteria  is  necessary  Co  insure  Chat  points  R  and  S  do 
not  fall  ouCside  points  A  and  B,  thus  maintaining  the  validity  of  Che 


inCerpolaCion.  From  Fig.  2-1 


which  becomes,’ 


< Qr”Qa> / < Qc’Qa <  XR~*A> / < *C~*A> 


Qr"Qc_^Qc~Qa) 


(2.23) 


(2.24) 


afCer  approximating  V+a  by  _+a  in  Eq.  (2.15).  This  approximation  is  valid 
where  V«a,  which  holds  for  sost  transient  problems  in  fluid  pipelines 
(13:58). 
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By  interpolations  similar  to  Eq.  (2.23) 


Qs“Qc-^Qc-Qb) 

(2.25) 

PR-Pc-r<Pc-PA> 

(2.26) 

pS“pC~f (Pc~pb) 

(2.27) 

where  f“aAt/Ax. 

To  calculate  the  final  solution  for  volume  flow  rate  and  pressure  at 
P  it  is  necessary  to  solve  six  equations.  The  six  are  Eqs.  (2.23) 
through  (2.27),  and  Eq.  (2.20)  and  (2.21).  The  addition  of  Eq.  (2.18) 
and  (2.19),  followed  by  the  solution  of  the  resulting  equation  for  the 
flow,  Qp,  gives: 

Qp-(CM+Cp)/2  (2.28) 

The  flow,  Qp,  may  then  be  substituted  into  either  Eq.  (2.18)  or  (2.19)  to 
obtain  the  pressure  at  that  point. 


III.  Boundary  Conditions 


The  equations  developed  in  the  lest  section  allow  the  calculation  of 
pressure  and  volume  flow  rate  at  any  point  P  within  the  pipe.  However, 
by  examining  Fig.  2-1,  the  end  points  of  the  grid  which  correspond  to  the 
pipe  under  consideration,  begin  influencing  the  interior  points  after  the 
first  time  step.  Therefore,  to  have  a  complete  solution  for  any  time 
after  the  first  time  step,  the  appropriate  boundary  conditions  become 
necessary. 

For  any  one  pipe,  only  one  characteristic  equation  is  available  at 
either  end,  as  seen  in  Fig.  3-1.  Referring  back  to  Fig.  2-1  and  the 
discussion  in  the  last  section,  Eq.  (2.19)  is  valid  at  the  upstream  boun¬ 
dary,  while  Eq.  (2.18)  holds  at  the  downstream  end.  With  only  the  one 
characteristic  equation  available  at  either  end,  a  second  equation  is 
required  to  complete  the  solution.  Thus,  to  find  the  pressure  and  the 
volume  flow  rate  at  the  end  of  a  pipe,  auxilary  equations  are  needed  to 
be  developed. 


Known  End  Conditions 

If  either  the  pressure  or  volume  flow  rate  at  an  end  of  the  pipe  is 
a  known  function  of  time,  F(t),  this  information  can  be  combined  with  the 
appropriate  characteristic  equation  to  fix  the  end  point  conditions.  For 


example,  if  the  pressure  at  Che  inlet  of  the  pipe  changes  in  a  known 
manner,  say  as  a  sine  wave,  the  boundary  condition  is 

F(  t)-Pp-P0+APsinwt  (3.1) 

where  u  is  the  frequency  and  AP  is  the  amplitude  of  the  wave.  With  the 
pressure  known  at  any  instant,  the  flow  at  the  inlet  is  then  determined 
by  direct  solution  of  Eq.  (2.19). 

Series  Connection 

From  Fig.  3-2,  recalling  Figs.  2-1  and  3-1,  it  can  be  inferred  Eq. 
(2.18)  is  available  for  pipe  1,  and  Eq.  (2.19)  is  available  for  pipe  2. 
However,  an  additional  equation  is  required  to  solve  each  characteristic 
equation.  The  continuity  expression  provides  one  equation,  while  a 
second  equation  is  obtained  by  equating  the  pressure  on  either  side  of 
the  junction  after  the  losses  are  assumed  to  be  negligible  (16:9.2). 

QP1"QP2  (3.2) 

Ppi*Pp2  (3.3) 

Along  with  Eq.  (2.18)  and  (2.19),  Eq.  (3.2)  and  (3.3)  provide  four 
equations  and  four  unknowns  which  may  then  be  solved  for  the  pressure  at 
the  junction. 

Pp-(Cp-CM)/(Bi+B2)  (3.4) 


The  remaining  unknowns  are  calculated  directly  from  the  appropriate 
equation. 


Fig.  3-2.  Series  Connection  Fig.  3-3.  Valve-in-line 

(16:9.3)  (16:9.5) 
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When  a  valve  or  orifice  is  located  between  two  pipes,  the  necessary 
equations  are  provided  in  much  the  same  manner  as  in  a  series  connection. 
Respectively,  Eq.  (2.18)  and  (2.19)  are  available  where  appropriate  as 
before,  and  the  continuity  expression,  Eq.  (3.2),  is  still  valid. 

However,  the  fourth  equation  is  supplied  by  the  valve  pressure-discharge 
characteristic,  r  (13:695).  For  positive  flow  (Fig.  3-3) 

T  »(Qp/Q0H  V<pPl-pP2>]*  (3*5) 

where  PQ  is  the  drop  in  pressure  across  the  fully  open  valve  at  the  ini¬ 
tial  flow  Q0,  and  APp  is  the  drop  in  pressure  across  the  partially  closed 
valve  at  an  instantaneous  flow  Qp.  By  solving  the  above  eqations 
simultaneously,  a  quadratic  in  Qp  results  which  may  be  solved  to  yield: 

Qp,-Cv(l/B1+l/B2)+(Cv2(l/B1+l/B2)2+2Cv(Cp/B1-CM/B2))i  (3.6) 

where  Cv=Q02r2/2P0.  If  there  is  flow  in  the  negative  direction  then  Eq. 
(3.5)  becomes 


JK 

fc 


I 


I 


r  *-(Qp/QoHP0/(pP2“pPl)l*  (3.7) 

which  yields  the  solution 

Qp*Cv(l/B1+l/B2)+[Cv2(l/B1+l/B2)2+2Cv(Cp/B1-CM/B2)]i  (3-8) 

Upon  examination  of  the  equations,  it  is  only  possible  to  have  negative 
flow  if  Cp/Bi-CM/B2  <  o.  Hence,  Eq.  (3.6)  is  valid  when  Cp/Bi-CM/B2  >  0, 
and  Eq.  (3.8)  is  used  when  Cp/Bi-Cm/B2  <  0.  With  the  flow  now  deter¬ 
mined,  Eq.  (2.18)  and  (2.19)  are  used  to  find  the  pressure,  Pp. 

Valve  at  Downstream  End  of  Pipe.  A  special  case  of  the  valve-in- 


line  is  a  valve 

ditions.  For  t 

i 

manner  (14:38) 

s 

r 

r“(Qp/Qo><po/Ppi>* 


(3.9) 
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Thus,  only  Eq.  (2.18),  Che  downstream  characteristic  equation,  is 
necessary  for  a  complete  solution.  This  simplification  results  in  Eq. 
(3.6)  reducing  to 

QP—  CV/B+[(CV/B)2+2CVCP]*  (3.10) 

While  the  corresponding  pressure  is  computed  from  either  Eq.  (2.18)  or 
(3.9). 

Cavitation 

During  transient  flows,  the  pressure  may  drop  below  the  vapor 
pressure  of  the  fluid.  When  this  occurs  the  fluid  undergoes  vaporiza¬ 
tion,  causing  vapor  bubbles  to  appear  in  the  flow.  The  phenomena  can  be 
treated  in  two  different  manners.  The  first  method  is  to  assume  the 
vapor  bubbles  are  dispersed  homogeneously  throughout  the  liquid.  The 
effective  wave  speed  resulting  from  this  two-phase  flow  is  then  treated 
as  a  function  of  pressure,  temperature,  and  vapor  volume  (14:136).  The 
alternate  method,  which  is  used  in  this  study,  lumps  the  free  gas  at  the 
computing  sections.  The  liquid  in  the  reaches  between  the  gas  volumes  is 
assumed  to  be  pure  liquid  without  free  gas.  This  assumption  allows  the 
use  of  a  constant  wavespeed  in  the  reaches  between  the  cavities  (17:49). 

When  the  pressure  falls  below  the  vapor  pressure,  the  cavity  boun¬ 
dary  condition  becomes: 

Pp»Pv  (3.11) 

where  Pv  is  the  vapor  pressure.  The  flow  upstream  of  the  cavity  is  then 
calculated  from  Eq.  (2.19) 

Qpy*Cp— Bp  (3.12) 

and  the  flow  downstream  by  Eq.  (2.18) 

Qp«CM+BPp  (3.13) 

12 


The  liquid  mess  conservation  is  maintained  by  applying  the  local 


continuity  relation  at  each  gas  volume. 


dVg/dt«Qd-Qu 


(3.14) 


Integration  of  Eq.  (3.14)  gives 


V “Vg+At 1 (Qp+Q)_(Qpu+Qu) 1/2  (3.15) 

where  Vg'  is  the  gas  volume  at  the  current  time  and  Vg  is  the  gas  volume 


from  the  previous  time  step  (14:137).  As  long  as  the  cavity  size  is 


positive,  vapor  pressure  persists.  Once  the  volume  becomes  less  than 


zero  the  gas  cavity  is  declared  to  have  collapsed  and  the  flow  calcula¬ 


tions  proceed  as  usual  for  an  interior  section  (17:49). 


However,  in  some  instances  the  system  to  be  analyzed  may  contain 


other  gases  dissolved  in  the  fluid.  This  phenomenon  is  common  in  rocket 


propulsion  where  fuel  and  oxidizer  tanks  are  often  pressurized  with  high 


pressure  helium  or  nitrogen.  Due  to  the  presence  of  this  other  gas  in 


the  fluid,  the  bubbles  which  occur  during  cavitation  will  contain  both 


the  vapor  of  the  fluid  and  the  gas  evolved  from  the  fluid.  Wylie  (17) 


developed  a  method  for  handling  the  additional  gas  release  in  the 


discrete  bubble  method.  The  eqations  in  terms  of  the  notation  of  this 


study  are: 


Pp»  -B+2(Pv+Pb)+[(Bx+2(Pv+Pb))2+8C4]* 


(3.16) 


where , 


Bi-[<^(CM-Cp)+Vg/2At  +  (  1-*)(0Q-0QU)  ]  /*B 


C4-Ci/2AtB<0 


C1-P*a0V 


In  which  a 0  is  the  void  fraction  at  some  reference  pressure,  P*,  while  V 


is  the  volume  of  the  adjacent  reach.  By  choosing  the  void  fraction  to  be 
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Accumulator 

In  this  study  an  accumulator  is  a  pressure  vessel  connected  to  the 
main  pipeline.  The  sealed  vessel  contains  a  layer  of  compressed  gas 
overlaying  the  liquid  inside.  The  conditions  which  exist  at  this  loca¬ 
tion  are  reminiscent  of  those  around  a  cavitation  bubble.  Thus,  Eq. 
(3.14)  is  once  again  used  to  express  the  continuity  relationship,  and  Eq. 
(2.18)  and  (2.19)  are  valid  where  appropriate.  However,  the  accumulator 
pressure  is  greater  than  vapor  pressure,  so  a  fourth  equation  is  needed 
to  complete  the  solution.  Dorsch  et  al .  (10:8)  obtained  the  necessary 
equation  by  introducing  the  compliance,  a',  of  the  gas  volume. 

a'-dVg/dP  (3.18) 

Making  use  of  the  chain  rule  from  calculus 

dVg/ dtm(dVg/dP) ( dP/dt )*a 1 dP/dt 


(3.19) 


uuuuuiww  wma 


I 

$ 

5 


•  •  < 


& 


Therefore 

dP/dt-(Qd-Qu)/a’  (3.20) 

Finally,  substituting  Eq.  (2.18)  and  (2.19)  into  Eq.  (3.20)  gives 

dP/dt-(CM+B2P-CP+B1P)/a'  (3.21) 

The  pressure  at  the  accumulator  is  now  a  known  function  of  time. 

Eq.  (3.21)  can  then  be  numerically  integrated  to  find  pressure  at  any 
time.  In  this  analysis  a  first  order  integration  is  used  to  give 

Pt+At^t*  At  ( CH4Vg+B2PaVg-Cpavg+BiPaVg)  /  a 1  (3.22) 

The  average  values  in  Eq.  (3.22)  are  determined  by  iteration  (see 
Appendix  A  for  a  listing  of  the  computer  code).  The  pressure,  Pp-Pt+Ac, 
so  Qp  above  and  below  the  accumulator  can  be  found  from  Eq.  (2.18)  and 
(2.19)  respectively. 


I 


V, 


3 

5 


V 


15 


I 

i 

In  this  section,  the  equations  developed  in  section  II  by  the  method  I 

( 

( 

of  characteristics  and  in  section  III  for  the  various  boundary  conditions  1 

I 

I 

were  used  to  model  several  pipe  systems  found  in  the  literature.  The  j 

systems  examined  began  vith  a  simple  one-pipe  system  and  moved  on  to  j 

systems  of  greater  complexity.  The  results  published  about  these  systems 
were  then  used  in  section  V  to  verify  the  accuracy  of  the  digital  computer 
program  written  for  this  study. 

Single  Pipe  Models 

As  mentioned  above,  the  first  systems  that  were  modeled  center 
around  a  single  pipe.  The  single  horizontal  pipe  was  combined  with 
various  upstream  and  downstream  boundary  conditions  to  verify  the 
accuracy  of  the  models  used  in  this  study  and  to  examine  the  effect  of 
such  boundary  conditions  on  a  system.  Note,  in  each  case  the  system  was 
assumed  to  have  achieved  steady  flow  before  any  changes  in  the  boundary 
conditions  were  imposed. 

Valve  Closure.  Wylie  and  Streeter  (14:39)  presented  the  simple  pipe 
flow  problem  of  a  sudden  valve  closure  at  the  end  of  a  single  pipe.  In 
this  case  the  system  consisted  of  a  single  pipe  with  a  reservoir  of  water 
at  the  upstream  end  and  a  valve  at  the  downstream  end.  The  flow  in  the 
pipe  was  steady  until  at  some  time,  t*0,  the  valve  began  to  close  such 
that  the  valve  pressure-discharge  characteristic,  Gq.  (3.9),  followed 
the  exponential  law 

T  -(l-t/tc)Eo  (4.1) 

where  tc  is  the  time  at  which  the  valve  is  fully  closed. 
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The  input  data  taken  from  Wylie  and  Streeter  for  this  system  were: 

D  *  0.5  m,  L  ■  600  m,  a  *  1200  m/s,  f  3  0.018,  PQ  3  1.4716  MPa,  Q0  3 
0.477  m^/s,  tc  *  2.1  s,  and  Em  =  1.5.  The  upstream  boundary  condition 
was  modeled  as  a  constant  pressure  reservoir.  With  the  pressure  at  the 
reservoir  a  known  function  of  time,  the  volume  flow  rate  was  determined 
by  Eq.  (3.2).  At  the  downstream  end  the  flow  was  discharging  through 
the  valve  to  ambient  conditions,  so  the  flow  through  the  valve  was 
calculated  by  Eq.  (3.11).  Hence,  the  pressure  at  the  downstream  end 
could  be  computed  with  either  Eq.  (3.9)  or  Eq.  (2.18). 

In  an  attempt  to  discover  the  effect  of  several  parameters  on  the 
severity  of  the  fluid  transients,  variations  were  introduced  into  the 
input  data.  The  variables  varied  appear  in  Table  I,  along  with  the  par- 


Table  I. 

Input  Variations  for 

Single  Pipe 

Closure 

i 

Run  No. 

tc  (s) 

Em 

f 

L  (m) 

1 

2.1 

2.0 

0.018 

600 

v\ 

•?: 

2 

2.1 

1.0 

0.018 

600 

3 

2.1 

0.5 

0.018 

600 

4 

1.5 

1.5 

0.018 

600 

> 

5 

2.5 

1.5 

0.018 

600 

6 

2.1 

1.5 

0.036 

600 

3 

;y 

7 

2.1 

1.5 

0.072 

600 

8 

2.1 

1.5 

0.0001 

600 

s'*, 

£ 

9 

2.1 

1.5 

0.018 

300 

10 

2.1 

1.5 

0.018 

150 

s 

11 

2.1 

1.5 

0.018 

75 

12 

2.1 

1.5 

0.018 

10 

■j 


1 


r 

w 


i 

,i 
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ticular  values  for  each  simulation.  The  remainder  of  the  variables  were 
kept  the  same  as  those  presented  above. 

Vaporous  Cavitation.  To  demonstrate  the  cavitation  model,  a  simple 
system  examined  by  Wylie  (17:50)  was  analyzed.  During  his  study  Wylie 
worked  with  two  systems  having  similar  configurations.  In  each  the 
pressure  at  both  ends  of  the  pipe  was  a  known  function  of  time.  The 
downstream  end  of  both  systems  was  connected  to  a  constant  pressure 
reservoir.  At  the  upstream  end  both  pipes  were  connected  to  a  pump, 
which  Wylie  represented  by  the  pressure  as  a  given  function  of  time. 
However,  the  two  systems  differed  as  to  the  physical  dimensions  of  the 
pipes,  as  well  as  to  the  properties  and  conditions  of  the  fluids  in  the 
systems . 

For  the  first  of  the  two  examples,  the  system  parameters  were: 

D  *  0.61  m,  L  -  3048  m,  f  =■  0.02,  Qp  *  0.89  m3/s,  a  *  981  m/s,  p  *  1000 
kg/m3,  and  Pv  *  -98.72  kPa  gage.  The  negative  pressure  resulted  from  the 
pressure  being  gage.  The  pressure  at  the  upstream  end  of  the  pipe  was 
dropped  linearly  from  an  initial  value  of  0.495  MPa  gage  to  zero  in  0.2 
8,  then  held  there  for  the  remainder  of  the  simulation  (17:51).  Thus, 
the  pressure  at  both  ends  of  the  pipe  was  known  for  any  time,  allowing 
the  volume  flow  rate  to  be  calculated  with  either  Eq .  (2.18)  or  (2.19). 

In  the  second  of  the  two  examples,  the  system  data  not  included  in 
Fig.  4-1  were:  Pm  1000  kg/m3,  and  Pv  ■  -99.11  kPa  gage.  As  in  the 
first  case,  the  downstream  boundary  condition  was  again  a  constant 
pressure  reservoir.  The  upstream  conditions  however  were  of  a  more 
complex  form  as  depicted  in  Fig.  4-1. 
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In  an  effort  to  further  verify  the  nodel  used  in  study,  two  addi¬ 
tional  cases  were  constructed  by  altering  the  original  input  data  given 
for  the  first  simulation.  The  new  hypothetical  systems  retained  the  same 
input  values  as  before  with  the  exception  of  the  vapor  pressure.  In  the 
one  case  the  vapor  pressure  was  set  at  ^--5.03  m,  while  in  the  other 
case  Hy«-100  m. 

Multipipe  Models 

Valve-in-line.  The  next  system  to  be  used  for  comparison  was  a 
system  analyzed  by  Swaf field  (13).  In  the  Swaf field  system  a  valve  was 
closing  in  the  middle  of  a  series  of  pipes  instead  of  at  the  end.  For 
that  investigation  Swaffield  used  aviation  kerosene  as  the  working  fluid 
not  water.  The  pipeline  itself  was  modeled  as  consisting  of  three  pipes 
of  the  equal  diameter  in  series,  where  both  the  upstream  and  the 
downstream  ends  of  the  system  were  connected  to  reservoirs.  The  valve  in 
question  was  situated  between  the  first  and  second  pipe,  and  as  in  the 
previous  valve  problem,  there  was  steady  flow  until  time,  t*0,  when  the 
valve  began  to  close.  The  parameters  which  describe  the  system  are  given 


in  Tables  II  and  III. 


Table  II.  Flow  Parameters  for  Valve-in-line 


Density, 

P *800  kg/m^ 

Wave 

Speed,  a*918  m/s 

tc  Downstream 

Pres.  Initial  flow 

Run 

(s) 

( kPa) 

(m ?/ s ) 

1 

0.081 

102 

0.00355 

2 

0.160 

222 

0.00355 

3 

0.140 

120 

0.00541 

Table  III.  Pipeline  Parameters  for  Valve-in-line 


Pipe 

Length 

(m) 

Diameter 

(m) 

f 

1 

5.80 

0.0508 

0.02 

2 

3.17 

0.0508 

0.02 

3 

6.70 

0.0508 

0.02 

Starting  with  the  upstream  end  of  the  pipeline,  there  was  once  again 
a  reservoir.  The  pressure  was  assumed  to  be  constant,  allowing  the 
volume  flow  rate  to  be  calculated  from  Eq.  (2.19).  Moving  downstream, 
the  next  boundary  condition  encountered  was  the  valve  between  pipes  1  and 
2.  Here,  as  shown  in  section  III,  the  flow  was  given  by  Eq.  (3.6),  where 
the  valve  pressure-discharge  characteristic  was  assumed  to  obey  Eq. 

(4.1).  The  pressure  on  the  upstream  and  downstream  face  of  the  valve 
could  then  be  calculated  from  Eq.  (2.18)  and  (2.19)  respectively.  At  the 
junction  of  pipes  2  and  3,  there  was  a  simple  series  connection.  The 
pressure  was  determined  by  Eq.  (3.4)  and  the  flow  was  then  calculated  once 
again  using  either  Eq.  (2.18)  or  (2.19).  The  final  boundary  condition 
was  the  constant  pressure  reservoir  at  the  downstream  end.  The  con¬ 
ditions  here  were  evaluated  by  substituting  the  known  pressure  once  more 
into  Eq.  (2.18). 


Valve  Closure.  A  second,  more  complex  downstream  valve  closure 
problem  was  also  presented  by  Wylie  and  Streeter  (14:60).  Like  the  first 
problem,  water  was  once  again  flowing  from  a  reservoir  at  the  upstream 
end  through  a  pipeline  to  a  valve  at  the  downstream  end.  In  this  example 
though,  the  pipeline  consisted  of  three  pipes  of  different  diameter  in 
series  between  the  reservoir  and  the  valve.  Furthermore,  the  valve 
pressure-discharge  characteristic,  Eq.  (3.9),  did  not  follow  an  exponen¬ 
tial  law,  instead,  r  follows  the  curve  in  Fig.  4-2. 

1 
r 

0.5 

0  I  2 

Time.  * 

Fig.  4-2.  Valve  characteristic,  r 

In  this  three  pipe  problem,  the  initial  flow,  Q0,  was  0.2  m^/s,  and 
the  initial  pressure  at  the  valve  equaled  981.4  kPa  gage.  The  parameters 
for  the  pipes  are  given  in  Table  IV. 


As  in  the  single  pipe  problem,  the  reservoir  was  assumed  to  be 
constant  pressure,  which  allowed  the  flow  to  be  calculated  at  that  point 
with  Eq.  (2.19).  Also,  the  flow  through  the  valve  was  again  calculated 
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using  Eq.  (3.10).  Therefore,  the  pressure  at  the  valve  could  be  calcu¬ 
lated  with  either  Eq.  (2.18)  or  (3.9).  However,  in  contrast  to  the  ori¬ 
ginal  single  pipe  closure  problem  there  were  now  two  interior  boundary 
conditions  as  well.  At  both  of  the  locations,  there  were  changes  in 
area.  The  pressure  at  such  junctions  is  calculated  by  Eq.  (3.4).  The 
flow  through  the  junction  could  then  be  determined  from  either  Eq.  (2.18) 
or  (2.19). 

Since  propellant  feedline  systems  are  rarely  a  single  pipe  or  a 
constant  diameter,  this  multipipe  valve  closure  problem  provided  an 
excellent  system  to  further  examine  the  effect  of  certain  parameters  on 
the  fluid  transients.  The  system  was  simple  enough  to  require  relatively 
short  run  times  on  the  computer,  about  10  CPU  seconds  on  the  Cyber.  Yet, 
at  the  same  time,  the  system  contained  features  that  are  present  in  the 
final  model  which  is  described  later. 

As  before  with  the  single  pipe  system,  while  keeping  the  other 
values  the  same  as  in  the  original  system,  the  variations  seen  in  Table  V 
were  introduced  into  the  input  data  one  at  a  time. 

Table  V.  Input  Variations  for  Series  Closure 


Run  No. 

Reservoir 
Pressure  (m) 

Closure 

1 

289 

original 

2 

400 

original 

3 

200 

original 

4 

289 

hard 

5 

289 

soft 

The  closure  designation  indicates  the  rapidity  of  the  closing  of  the 
valve.  The  original  closure  was  the  valve  characteristic  in  Fig.  4-2. 


During  a  hard  closure,  the  value  of  r  decreased  from  1.0  to  0.1  in  the 
first  0.6  s,  then  went  to  zero  in  the  next  0.6  s.  And  lastly,  in  the 
soft  closure,  r  dropped  from  1.0  to  0.4  during  the  first  0.6  s,  went  from 
0.4  to  0.15  in  the  next  0.6  s,  then  fell  to  zero  during  the  final  0.6  s. 
Thus,  the  characteristic  curve  for  the  hard  closure  not  only  had  a 
steeper  slope,  but  also  a  cut-off  time  that  was  0.6  s  shorter  then  the 
other  two  curves. 

In  addition  to  running  each  case  with  water  as  the  working  fluid, 
runs  were  also  made  with  rocket  fuels  RP-1,  a  fuel  very  similar  to  the 
aviation  kerosene  used  by  Swaf field,  and  liquid  hydrogen.  However,  in 
order  to  have  a  complete  set  of  input  data  for  these  addition  fluids, 
several  assumptions  were  made.  The  first  of  the  assumptions  was  asso¬ 
ciated  with  the  wave  speeds.  When  a  wave  speed  is  not  explicitly  stated 
in  the  input  list,  the  program  calculates  the  wave  speed  as  follows 
(14:58) 

a2*K/p( 1+DK+Ee)  (4.2) 

where  the  bulk  modulus,  K  and  the  density,  p,  are  fluid  properties,  while 
the  Young' 8  modulus,  E,  the  pipe  dameter,  D,  and  the  pipe  thickness,  e, 
are  properties  of  the  pipe. 

The  wave  speed  for  water  was  given,  along  with  the  diameter  of  the 
pipes  for  the  series  closure  by  Wylie  and  Streeter  (14:60),  but  nothing 
else.  Therefore,  the  water  was  assumed  to  have  a  density  of  1000  kg/m^ 
as  in  the  single  pipe  cavitation  problems.  Given  this  density,  a  bulk 
modulus  of  2.016x10^  N/m2  was  obtained  from  Appendix  B  of  Dehoff  (2). 
Also,  the  pipes  were  assumed  to  have  a  Young's  modulus  of  7.24xl010  N/m2, 
the  value  reported  by  Swaffield  (13:698)  for  his  system.  With  everything 
but  one  variable  now  known  in  Eq.  (4.2),  it  was  possible  to  determine  a 


thickness  for  each  pipe.  Substituting  the  appropriate  values  into  the 
equation  gave  ei*0.021  m,  e2M0.014  m,  and  e3“0.01  m,  where  ex  is  the 
thickness  of  pipe  x.  Thus,  a  complete  description  of  the  pipe  system  was 
produced.  As  for  the  fluid  properties,  they  are  listed  below: 


RP-1 

lh2 

Temperature  (°F) 

89.10 

-423.2 

Density  (kg/m-*) 

800.00 

72.1 

Bulk  Modulus  (10®  Pa) 

12.34 

0.1 

Vapor  Pressure  (Pa) 

700.00 

127,000 

The  second  set  of  assumptions  centered  around  the  friction  factor. 

In  equation  form  the  friction  factor  is  given  by 

f*2APfD/pLV2  (4.3) 

where  APf  *  pressure  loss  due  to  friction,  V  *  mean  flow  velocity,  and 
L*pipe  length.  For  the  simulations  where  the  RP-1  and  LH2  replaced  water 
everthing  in  Eq.  (4.3)  but  the  density  remained  constant.  Thus,  a  fric¬ 
tion  factor  for  the  other  two  fluids  was  obtained  by  multiplying  the 
friction  factor  of  water  by  a  ratio  of  the  density  of  the  water  to  the 
density  of  the  other  fluid.  However,  this  process  produced  a  very  high 
friction  factor  for  LH2,  so  a  second  run  with  LH2  as  the  fluid  was 
included  for  each  case  of  the  simulation.  In  this  second  run  LH2  was 
given  the  same  friction  factor  as  water  to  determine  how  strongly  the 
friction  was  affectiug  the  transients  in  the  LH2. 

Saturn  V  Feedline.  The  final  system  analyzed  was  chosen  in  order  to 
test  the  routines  developed  in  this  investigation  on  a  real  aerospace 
system.  The  system  chosen  was  one  that  was  reported  on  by  Brod  (18). 
During  the  Apollo  program,  Boeing  constructed  a  test  installation  to 
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simulate  the  LOX  suction  duct  of  the  Saturn  V’s  F-l  engine.  In  the 
course  of  the  testing  performed  wich  this  apparatus,  some  experiments 
were  conducted  which  consisted  of  closing  the  main  drain  valve  and 
observing  the  fluid  transients  (18:15).  Unfortunately,  the  results  from 
these  experiments  were  no  longer  available  from  Boeing.  Thus,  the 
results  of  this  study  for  this  simulation  were  of  a  more  qualitative 


nature , 


Table  VI.  Feedline  Configuration  Data 
Pipe _ Length  (m)  Dia.  (m)  Wave  Speed  (m/s) 


0.508 

0.508 


0.508 


0.457 


0.500 


0.457 

0.457 


0.457 


0.457 


0.462 


0.640 


0.462  ' 


0.457 


0.457 


0.457 


770.7 


770.7 


770.7 


862.1 


658.4 


862.1 


1165.3 


1165.3 


862.1 


207.8 


177.5 


207.8 


800.2 


1330.0 


1330.0 


In  general,  this  problem  was  related  to  the  multipipe  example  taken 
from  Wylie  and  Streeter  (14).  In  both,  there  was  a  reservoir,  a  series 
of  pipes,  and  a  valve  at  the  downstream  end.  In  addition,  both  were 
modeled  as  horizontal  pipelines  for  this  study,  however,  the  similarities 
end  there.  The  F-l  feedline  was  comprised  of  15  pipes  in  all.  The  data 
for  the  pipes  was  given  in  Table  VI.  Furthermore,  in  addition  to  more 
pipes  in  the  system,  there  were  a  few  additional  interior  boundary  con¬ 
ditions.  In  the  junction  between  pipes  7  and  8  was  a  device  called  the 
prevalve.  Initially  the  prevalve  was  a  simple  visor  valve.  But,  as 
part  of  the  experiments  at  the  test  installation,  the  valve  cavity  was 
evacuated  and  pressurized  with  helium.  In  this  way  the  engineers  hoped 
to  control  the  fluid  transients  by  essentially  creating  an  accumulator  at 
that  point.  Therefore,  for  the  purposes  of  this  analysis  the  prevalve 
was  modeled  as  an  accumulator  when  there  was  a  volume  of  gas  present. 

The  second  new  boundary  condition  was  the  pulser  valve  between  pipes 
13  and  14.  To  model  this  component  an  additional  boundary  condition  was 
developed.  Recalling  similar  circumstances  in  Figs.  3-2  and  3-3,  it  was 
apparent  from  Fig.  4-3  that  Eq.  (2.18)  and  (2.19)  were  again  available. 
The  continuity  equation  for  this  junction  is  of  the  form 

Qpu*Q+Qp  (4.4) 

where  the  flow  out  of  the  pulser  valve  was  assumed  to  be  of  the  form 

Q“Qv|sinwt|  (4.5) 

The  final  relationship  required  to  form  a  closed  solution  was  taken  from 
a  work  by  Thorley  (16:9.3).  In  that  paper,  he  proposed  that  in  a 
branching  connection  such  as  this  a  common  pressure  could  be  assigned  at 
the  junction,  as  long  as  any  losses  were  neglected  as  minor.  Thus,  Eqs . 
(2.18),  (2.19),  and  (4.5)  were  substituted  into  the  continuity  equation 
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and  solved  for  Che  pressure  Co  give 

Pp“(Cp-Cjr"Qv j s inwC |  )/(Bj+B2)  (4.6) 

The  pressure  was  Chen  subsciCuCed  back  inCo  Eq.  (2.18)  and  (2.19)  Co  find 
che  flow  aC  those  secCions. 

The  remaining  incerior  boundary  locaCed  between  pipes  14  and  15  also 
had  a  boundary  condition  other  than  a  series  connection.  In  order  to 
complete  Che  simulation  of  the  LOX  suction  duct,  a  second  accumulator  was 
included  to  approximate  cavitation  at  the  inlet  of  the  turbopump. 

In  this  system,  the  parameters  reported  by  Dehoff  (2)  and  Brod  (18) 
were  as  follows:  Q0*1.7  m^/s,  P0  at  the  valve*621  and  896  kPa,  p*1000 
kg/m^,  w*0  to  27  Hz,  Qv*0.071  m^/s,  Vg  of  the  prevalve  accumulator*© , 

30,  60,  and  90  percent  of  the  total  accumulator  volume,  valve  cut-off 
times*1.25  to  4.0  a,  and  Vg*0.0129  m^  at  the  pump  inlet  when  cavitation 
was  included.  The  specific  values  arbitrarily  chosen  for  each  parameter 
which  had  a  range  of  values  is  shown  in  Table  VII. 
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Table  VII.  Input  Parameters  for  F-l  Feedline 


Run  no. 

Valve(kPa) 

tc  (s) 

Em 

(Hz) 

Vol . (m^) 

cav 

1 

621 

1.5 

1.2 

0.0 

0.000 

no 

2 

621 

3.0 

1.2 

0.0 

0.000 

no 

3 

621 

1.5 

1.0 

0.0 

0.000 

no 

4 

621 

1.5 

1.5 

0.0 

0.000 

no 

5 

896 

1.5 

1.2 

0.0 

0.000 

no 

6 

621 

1.5 

1.2 

2.5 

0.000 

no 

7 

621 

1.5 

1.2 

10.0 

0.000 

no 

8 

621 

1.5 

1.2 

25.0 

0.000 

no 

9 

621 

1.5 

1.0 

2.5 

0.000 

no 

10 

896 

1.5 

1.2 

2.5 

0.000 

no 

LI 

896 

1.5 

1.0 

10.0 

0.000 

no 

12 

621 

1.5 

1.0 

0.0 

0.000 

yes 

13 

621 

1.5 

1.0 

0.0 

0.018 

yes 

14 

621 

1.5 

1.0 

0.0 

0.036 

yes 

15 

621 

1.5 

1.0 

0.0 

0.072 

yes 

16 

896 

1.5 

1.0 

0.0 

0.000 

yes 

17 

896 

1.5 

1.0 

0.0 

0.054 

yes 

18 

621 

1.5 

1.0 

2.5 

0.054 

yes 

V .  Results 


Single  Pipe  Verification  for  Valve  Closure 

Some  of  the  published  results  used  for  verification  were  in  terms  of 
piezometric  head  instead  of  pressure*  Since  the  program  uses  equations  in 
terms  of  pressure,  a  conversion  was  required.  To  obtain  the  conversion 
the  equation,  P“Pg(H-z),  was  solved  for  H.  Thus,  when  necessary,  the 
resulting  conversion  was  incorporated  into  the  output  statement  of  the 
program,  making  use  of  FORTRAN'S  ability  to  perform  mathematical  opera¬ 
tions  on  the  variable  to  be  output. 

Simple  Valve  Closure.  The  first  system  to  be  analyzed  was  the 
single  pipe  valve  closure  presented  by  Wylie  and  Streeter  (14).  The 


results  published  by  Wylie  and  Streeter  for  this  problem  appear  in  Fig- 


5-1.  The  agreement  between  the  published  results  and  the  results  pre¬ 


dicted  by  the  program,  which  appear  in  Fig.  5-2,  was  very  close.  For  the 


pressure,  both  results  depicted  the  piezometric  head  increasing  to  a  peak 


of  about  285  m  at  the  valve  as  the  decreasing  area  generated  compression 


waves  which  traveled  upstream.  The  compression  waves  reflected  from  the 


constant  pressure  boundary  at  the  reservoir,  becoming  expansion  waves 


which  propagated  back  toward  the  valve.  Since  vapor  pressure  was  not 


reached,  the  expansion  waves  which  returned  to  the  valve  in  2L/aI>1.0  s 


were  stronger  than  the  compression  waves  being  generated  at  the  valve, 


dropping  the  head  pressure  at  that  point.  Here,  the  strength  of  a  wave 


was  equated  with  the  pressure  difference  across  a  wave.  Once  the  valve 


was  fully  closed,  the  waves  in  the  pipe  simply  reflected  between  the 


free  boundary  at  the  reservoir  and  a  wall  at  the  valve  with  a  period  of 


4L/a«2.0  s. 


The  predictions  of  the  pressure  matched  very  closely  when  compared 
point  by  point.  So,  in  order  to  further  explore  the  effect  of  certain 
system  parameters  on  the  resulting  fluid  transients,  the  variations 
listed  in  Table  I  were  introduced  into  the  system.  The  plots  of  the 
results  predicted  for  the  various  changes  appear  in  Figs.  5-3  through 
5-6. 

Effects  of  Closure  Curve.  The  first  of  the  parameters  examined  was 
the  valve  closure  characteristic,  r.  In  Fig.  5-3,  appear  the  results 
produced  by  varying  the  exponent,  Eq,  from  Eq.  (4.1),  the  equation  used 
to  represent  the  valve  closure.  Not  only  did  the  magnitude  of  the  first 
peak  change  with  changes  in  Eq,  but  the  slope  of  the  curve  as  the  peak 
was  approached  changed  with  Eq  as  well,  this  was  as  expected  though.  The 
timing  of  the  reflections  should  have  been  and  was  the  same  as  seen  in 
Fig.  5-2.  The  change  during  these  runs  was  expected  in  the  strength  of 
the  waves  at  a  particular  time  during  the  closure.  It  was  anticipated 
that  when  the  closure  was  initially  rapid  then  slowed,  that  the  pressure 
would  peak  more  quickly  than  when  the  closure  was  initially  slow  then 
became  increasingly  quicker.  Further,  it  was  anticipated  that  the  peak 
magnitude  would  increase  as  Eq  became  larger  or  smaller,  since  in  either 
case  the  closure  was  approaching  an  instantaneous  cut-off. 

Starting  with  an  Eq  greater  than  one,  the  valve  closed  rapidly  at 
first  but  slowed  as  the  cut-off  time  approached.  Consequently,  the 
strongest  compression  waves  were  generated  during  the  start  of  the  clo¬ 
sure,  while  the  weakest  occurred  at  the  end.  This  meant  that  when  the 
waves  returned  from  the  reservoir  at  2L/a*1.0  s,  the  reflected  expansion 
waves  were  stronger  than  the  compression  waves  currently  being  generated 
by  the  closing  valve.  So,  the  pressure  began  to  drop.  At  4L/as2.0  s  the 
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compression  waves  which  resulted  from  a  second  reflection  at  the  reser¬ 
voir  began  to  arrive.  This  third  set  of  waves  combined  with  the  wave 
trains  already  at  the  valve  to  decrease  the  rate  at  which  the  pressure 
decreased.  Then,  at  t»2.1  s,  the  valve  was  totally  closed  and  all  that 
remained  was  the  interaction  of  the  reflecting  wave  trains  as  they  travel 
back  and  forth. 

When  Eq  was  one,  the  valve  closed  at  a  constant  rate,  generating 
compression  waves  which  diminished  slightly  in  strength  as  the  flow  went 
to  zero.  Therefore,  when  the  expansion  waves  began  to  arrive  after  the 
first  reflection,  there  was  essentially  a  cancellation  between  them  and 
the  compression  waves  being  generated  at  that  time.  This  cancellation 
resulted  in  the  relatively  flat  portion  of  the  curve  observed  between  1.0 
and  2.0  s.  Then  at  t*2.0  s  the  second  reflection  from  the  reservoir 
arrived,  reducing  tht  rate  of  the  pressure  drop.  But  at  t*2.1  s,  the 
valve  was  completely  closed,  thus  compression  waves  were  no  longer  being 
generated.  And  so,  the  expansion  waves,  which  had  not  been  as  greatly 
attenuated  by  friction  as  the  compression  waves,  produced  a  pressure  drop 
that  continued  until  t»3.1  s,  when  the  last  of  the  wave  train  from  the 
first  reflection  finally  passed.  Then  once  more  all  that  remained  was 
the  wave  train  reflecting  back  and  forth. 

Finally,  for  Eq  less  than  one,  the  valve  shut  slowly  at  first  then 
closed  with  increasing  rapidity  as  the  cut-off  time  approached. 
Consequently,  the  strength  of  the  compression  waves  generated  at  the 
valve  increased  with  time.  Thus,  when  the  wavefront  returned  from  the 
first  reflection  at  1.0  s,  the  expansion  waves  were  not  as  strong  as  the 
compression  waves,  and  the  pressure  continued  to  rise,  but  at  a  slower 
rate.  However,  at  2.0  s  the  compression  waves  from  the  second  reflection 


arrived,  further  increasing  the  rate  of  the  pressure  rise.  But,  once 
more  the  valve  was  totally  closed  at  2.1  s,  allowing  the  expansion  waves 
from  the  first  reflection  to  dominate,  resulting  in  a  steep  pressure 
drop.  And  once  more,  all  that  remained  after  this  time  was  the  wave 
train  repeatedly  reflecting  within  the  pipe. 

Effects  of  Valve  Cut-off  Time.  Moving  now  to  the  second  parameter, 
the  length  of  time  required  for  the  valve  to  totally  close  was  examined 
for  its  effect  on  the  transient.  In  Fig.  5-4,  the  result  of  the  origi¬ 
nal  simulation  with  a  cut-off  time  of  2.1  s  was  compared  with  the  results 
of  a  simulation  where  the  cut-off  time  had  been  reduced  to  1.5  s,  as  well 
as  a  simulation  where  the  cut-off  time  was  increased  to  2.5  s.  For  the 
case  of  the  shorter  cut-off  time,  the  valve  closed  more  suddenly  than  in 
the  original  case.  This  more  sudden  closure  generated  stronger  compres¬ 
sion  waves,  hence  produced  a  steeper  pressure  rise  during  the  first 
second,  as  seen  in  Fig.  5-4.  Since  was  greater  than  one,  at  2L/a»1.0 
s,  the  expansion  waves  from  the  first  reflection  arrived  cauing  a  drop  in 
pressure.  However,  in  this  case  the  valve  was  fully  closed  by  1.5  s,  so, 
the  expansion  waves  did  not  encounter  any  other  waves  until  t*2.0  s  when 
the  compression  waves  from  the  second  reflection  of  the  wave  train 
arrived.  This  period  of  pure  expansion  dropped  the  pressure  near  vapor 
pressure  before  the  arrival  of  the  second  reflection  could  reverse  the 
trend.  But  ultimately,  this  case  settled  into  the  same  pattern  as  before 
with  pressure  oscillations  occuring  with  a  period  of  4L/aa2.0  s. 

In  the  case  of  a  cut-off  time  longer  than  tca2.1  s,  the  closure  was 
more  gradual  than  the  original,  generating  weaker  compression  waves  as  a 
result.  As  seen  in  Fig.  5-4,  these  weaker  waves  produced  a  slower 
pressure  rise  which  lasted  until  the  expansion  waves  from  the  first 


reflection  arrived  at  1.0  s.  When  the  compression  waves  from  the 
second  reflection  arrived  at  2.0  s  the  valve  was  still  closing,  thus  the 
two  sets  of  weaker  compression  waves  combined  to  reduce  the  rate  of  the 
pressure  decrease  produced  by  the  expansion  waves.  Then  at  2.5  s  the 
valve  was  totally  closed.  But,  even  though  compression  waves  were  no 
longer  generated  at  the  valve,  the  pressure  began  to  rise  due  to  the 
difference  in  wave  strengths  caused  by  the  relative  difference  in  times 
during  the  closure  at  which  the  waves  currently  at  the  valve  were 
generated.  Then,  at  3.0  s  expansion  waves  from  the  third  reflection 
began  to  arrive,  decreasing  the  rate  at  which  the  pressure  increased 
towards  the  second  peak.  With  the  valve  closed,  all  that  remained  after 
this  time  were  the  waves  reflecting  between  the  valve  and  the  reservoir. 

Effects  of  Friction  Factor.  The  next  parameter  investigated  for  its 
effect  on  the  pressure  at  the  valve  during  the  fluid  transient  was  the 
friction  factor  of  the  pipe.  The  results  of  changing  the  friction  factor 
are  plotted  in  Fig.  5-5.  As  the  friction  factor  increased,  the 
attenuation  due  to  the  friction  increased.  This  attenuation  was  apparent 
in  several  ways.  First,  the  initial  pressure  at  the  valve  decreased  with 
increasing  values  of  the  friction  factor.  Second,  the  rate  at  which  the 
pressure  oscillations  damped  out  increased  as  the  friction  factor 
increased.  And  third,  despite  an  increasing  difference  between  the  ini¬ 
tial  pressure  and  the  peak  pressure,  the  peak  magnitude  of  the  transient 
decreased  as  the  friction  factor  increased. 

Structurally,  the  peak  magnitude  is  the  most  important  factor. 
However,  in  this  study  the  initial  valve  pressure  was  of  more  interest 
due  to  its  relationship  with  the  timing  of  the  peak  magnitude  as  well  as 
the  timing  and  magnitude  of  the  secondary  peaks.  When  the  friction  fac- 


tor  was  small  there  was  little  attenuation  of  the  flow,  so  the  pressure 
gradient  between  the  reservoir  and  the  valve  was  very  slight.  Thus,  once 
the  valve  began  to  close,  the  pressure  at  the  valve  rose  until  the  first 
reflection  arrived  at  1.0  s.  At  this  time  the  returning  expansion  waves 
began  to  reduce  the  volume  flow  rate  almost  immediately  since  the 
pressure  gradient  was  small.  Afterwards,  the  wave  interactions  were 
like  those  described  in  the  original  case,  with  waves  reflecting  back 
and  forth  between  the  valve  and  the  reservoir. 

But,  as  the  friction  factor  and  the  attenuation  increased,  the  ini¬ 
tial  pressure  at  the  valve  decreased,  creating  a  steeper  pressure  gra¬ 
dient.  As  the  gradient  steepened,  the  strength  of  the  compression  at  the 
valve  increased.  Hence,  the  pressure  at  the  valve  rose  faster  initially. 
Also,  once  the  reflection  returned,  the  expansion  waves  took  longer  to 
reverse  the  flow,  causing  what  appears  as  a  displacement  in  time  of  the 
pressure  peaks.  For  the  case  of  £*0.072,  the  gradient  was  such  that 
despite  the  arrival  of  the  expansion  waves,  the  pressure  continued  to 
rise  significantly  before  the  flow  was  finally  reversed  and  the  pressure 
began  to  drop.  Further,  since  the  compression  strengthened  as  the  fric¬ 
tion  factor  increased,  the  subsequent  reflections  of  the  wave  train  were 
stronger  as  well.  This  meant  that  even  though  the  peak  magnitude 
decreased  with  increasing  friction  factor,  the  increasing  AP  caused  the 
the  magnitude  of  the  secondary  peaks  to  be  greater  initially,  but  damp 
out  at  a  higher  rate  due  to  the  increased  friction. 

Effects  of  Pipe  Length.  The  final  parameter  examined  for  its  effect 
on  the  transient  was  the  actual  length  of  the  pipe  under  consideration. 
The  results  produced  by  altering  the  length  of  the  pipe  appear  in  Fig. 
5-6.  As  the  pipe  becomes  shorter,  the  time  required  for  the  waves  to 
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traverse  the  length  of  the  pipe  decreases.  Consequently  when  L=300  m, 
the  first  reflection  arrived  at  the  valve  in  0.5  s,  while  the  system  had 
a  period  of  1.0  s.  For  La150  m,  the  reflection  arrived  in  0.25  s,  while 
the  system  period  was  0.5  s,  and  so  on  as  the  pipe  length  decreased. 
Moreover,  since  the  expansion  waves  returned  to  the  valve  with  increasing 
quickness,  the  pressure  had  increasingly  less  time  to  be  influenced  by 
the  compression  due  to  the  closure.  Thus,  the  magnitude  of  the  peak 
pressure  decreased  as  the  pipe  length  decreased. 

Aside  from  the  effects  already  mentioned,  the  reduced  travel  time 
had  two  other  effects.  First,  the  relative  strength  of  any  two  waves  was 
dependent  on  the  points  in  time  during  the  closure  at  which  each  of  the 
waves  was  generated,  since  the  closure  characteristic  was  nonlinear. 
Again,  strength  was  equated  to  the  pressure  difference  across  a  wave. 
Therefore,  the  shorter  the  travel  time,  the  closer  the  compression  waves 
currently  generated  at  the  valve  were  in  strength  to  the  reflected  expan¬ 
sion  waves  returning  from  the  reservoir.  Second,  as  the  travel  time 
decreased,  the  time  interval  between  reflections  decreased  as  well.  This 
meant  the  number  of  reflections  occuring  before  the  valve  was  totally 
closed  increased. 

The  significance  of  these  two  effects  was  that  the  combination  of 
the  two  produced  the  differing  responses  observed  in  the  time  interval 
between  the  arrival  of  the  first  reflection  and  the  complete  closure  of 
the  valve.  When  the  pipe  length  was  such  that  2L/a  was  a  significant 
percentage  of  the  closure  time,  only  a  few  reflections  occurred  during 
closure  and  the  relative  strengths  of  the  current  and  the  reflected  waves 
were  very  different.  Hence,  as  was  seen  in  Fig.  5-6,  for  L*300  m  the 
pressure  dropped  rather  rapidly  once  the  reflection  arrived  at  0.5  s.  In 
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addition,  it  via  easy  to  distinguish  the  additional  reflections  every  0.3 
s  after  that  point.  However,  as  the  pipe  became  shorter,  more  reflec¬ 
tions  of  about  the  same  relative  strength  began  interacting  at  the  valve 
more  quickly.  The  result  was  essentially  a  cancellation  among  all  of  the 
reflections,  which  left  the  relative  difference  between  the  strengths  of 
the  primary  compression  and  the  first  reflection  as  the  main  driving 
force  of  the  pressure  changes  observed  during  the  valve  closure. 


Single  Pipe  Verification  for  Cavitation 

First  Simulation.  In  the  first  of  the  cavitation  simulations  the 


upstream  pressure  was  dropped  from  495  kPa  to  zero  and  held  there.  The 
results  published  by  Wylie  for  this  system  are  shown  in  Fig.  3-7,  while 
Fig.  5-8  and  5-9  show  the  results  of  the  program  developed  in  this  study. 
In  general,  all  the  curves  agreed  with  the  expected  results.  That  was, 


the  pressure  downstream  of  the  "pump"  dropped  to  vapor  pressure  as  the 


Fig.  5-7. 


Wylie's  results  for  two  locations  (L*3048  m,  a*981  m/ s ,  D*0.61 
Qo-0.89  s,  Hy— 10.06  m) 


expansion  waves  generated  by  the  change  in  upstream  conditions  traveled 
down  the  pipe.  Further,  as  was  expected,  the  reflection  of  those  waves 
returned  at  some  time  greater  than  2L/a“4.56  s  as  indicated  by  the 
location  of  the  major  pressure  surge.  The  surge  itself  was  due  to  the 
expansion  waves  reflecting  as  compression  waves  from  the  constant 
pressure  boundary  at  the  reservoir,  then  propagating  back  upstream.  The 
delay  in  the  return  of  the  waves  seen  in  the  results  was  introduced  by 
the  formation  of  vapor  within  the  pipe  system.  As  discussed  by  Wylie 
(14,  17),  the  effective  wave  speed  is  inversely  related  to  the  mass  of 
the  vapor  present.  Thus,  as  vapor  was  released  the  effective  wave  speed 
decreased,  creating  the  delay  in  the  occurrence  of  the  pressure  surge  at 
the  station  indicated. 

Figs.  5-7,  5-8,  and  5-9  agreed  very  closely  as  to  the  time  and 
approximate  magnitude  of  the  major  events.  Unexpectedly,  the  exact 
magnitude  of  the  program  results  was  obscured  by  numerous  pressure  spikes 
in  the  solution.  The  period  of  the  spikes  was  related  to  the  distance 
between  the  downstream  end  of  the  cavitation  bubble  and  the  end  of  the 
pipe.  The  cause  of  those  spikes  however,  will  be  discussed  later. 

Second  Simulation.  For  the  second  of  the  cavitation  simulations, 
the  upstream  conditions  became  those  seen  in  the  bottom  curve  in  Fig. 
5-10.  The  results  Wylie  published  for  this  slightly  more  complex 
interaction  appear  in  Fig.  5-10  a3  the  upper  two  curves.  The  results  of 
this  study  follow  in  Fig.  5-11  and  5-12.  In  this  case  the  change  in  the 
upstream  boundary  condition  generated  both  expansion  and  compression 
waves.  The  comparison  of  the  curves  was  quite  good,  considering  the 
approximation  of  the  upstream  pressure  that  was  used  for  the  program 
results.  Both  sets  of  results  agreed  that  the  first  peak  was  delayed  due 
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Co  Che  pressure  falling  Co  vapor  pressure.  Also,  boCh  sees  indicated  Che 
syscem  had  a  period  of  2.25  s  once  Che  cavities  had  collapsed,  a  period 
which  was  half  of  the  expected  4L/a»4.5  s.  This  anomaly  was  explainable 
by  realizing  both  Che  expansion  and  the  compression  which  occurred  ac  Che 
upscream  end  were  of  approximately  the  same  magnitude  and  were  seperated 
in  cime  by  almost  exactly  L/a*l.l  s.  Thus,  Che  expansion  and  Che 
compression  seperaCed  by  this  cime  interval  caused  the  system  to  act  as 
if  the  pipe  was  only  half  of  its  actual  length.  But,  even  chough  Che 
magnitude  of  Che  peaks  agreed  well,  there  were  scill  minor  discrepancies 
in  Che  exact  shape  and  magnitude. 

Further  Verification.  To  verify  that  the  discrete  vapor  model  was 
indeed  working  correctly  in  light  of  the  unrealistic  pressure  spikes,  the 
first  cavitation  simulation  was  repeated,  but  with  two  different  vapor 
pressures.  Theoretically,  if  two  fluids  with  different  vapor  pressures 
were  exposed  to  the  same  drop  in  pressure,  the  fluid  with  the  higher 
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Fig.  5-12.  Program  result*  for  simulation  2  (x-580  m,  L*1450  m,  D-0.1  m, 
a-1290  m/s, Qo-0. 0158  m3/s,  Hv— 10.1  m) 


l 


v 


43 


r 


vapor  pressure  would  release  a  greater  amount  of  vapor.  Further,  the 
larger  mass  of  vapor  would  produce  a  larger  decrease  in  the  effective 
wave  speed,  given  the  same  initial  wave  speed  in  both  fluids.  A  vapor 
pressure  of  -100  m  was  first  chosen  to  insure  that  the  fluid  remained  a 
liquid  throughout  the  simulation.  This  choice  of  vapor  pressure  made  it 
possible  to  verify  that  the  model  was  correctly  handling  the  wave  propa¬ 
gation  and  attenuation.  For  example,  the  pressure  surge  at  x*813  m 
downstream  from  the  inlet  should  occur  2L/a=4.56  s  after  the  first  expan¬ 
sion  wave  passes  that  point.  That  is  the  time  required  for  a  wave  tra¬ 
veling  981  m/s  to  propagate  downstream  and  return.  For  the  results  in 
Fig.  5-13,  the  initial  wave  passed  x*813  m  at  about  0.8  s,  the  surge 
began  at  5.4  s.  This  made  a  difference  of  4.6  s,  which  is  within  one 
percent  of  the  predicted  value.  As  for  the  attenuation  of  the  waves,  the 
decrease  in  the  strength  of  the  waves  could  be  seen  by  the  manner  in 
which  the  curves  converged  towards  zero. 

For  the  case  of  a  vapor  pressure  higher  than  the  original  case,  the 
previous  argument  was  still  true.  That  is  to  say,  the  higher  the  vapor 
pressure,  the  greater  the  amount  of  vapor  released,  and  the  lower  the 
wave  speed  must  be  for  a  given  pressure  drop.  Thus,  when  the  vapor 
pressure  was  raised  to  -5.03  m  gage  in  the  cavitation  simulation,  the 
pressure  surge  was  expected  to  occur  at  some  later  time.  A  comparison  of 
Fig.  5-8,  the  results  for  ^*-10.06  m,  and  Fig.  5-14,  the  results  for 
Hv*-5.03  m,  indicated  that  the  pressure  surge  was  indeed  delayed, 
verifying  that  much  of  the  model.  However,  the  results  were  again 
plagued  by  the  nonphysical  pressure  fluctuations. 

Pressure  Spikes.  These  fluctuations  and  spikes  were  attributed  to 
the  method  by  which  the  vapor  was  accounted  for  in  the  flow.  By  lumping 


Fig.  5-13. 


Program  results  for  simulation  1  (L*3048  m,  a*981  m/s,  D“0.61  m 
Qo*0 .89  s,  Hy*-100  m) 


Fig.  5-14.  Program  results  for  simulation  1  (**813  m,  L“3048  m,  a*981  m/s, 
D*0.61  m,  Qo*0.89  s,  Hy— 5.03  m) 
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Che  gas  bubbles  at  Che  compuCing  secCions,  numerical  oscillacions  in  Che 
magnicude  of  Che  pressure  occurred.  These  nonphysical  pressure  oscilla¬ 
cions  arose  for  several  reasons.  FirsC  when  Che  gas  volume  was  lumped  ac 
a  compuCing  secCion,  ic  effecCively  produced  an  inCernal  conscanc 
pressure  boundary.  Hence,  when  mulciple  cavicies  occurred,  chere  were 
mulciple  reflecCions  of  Che  waves  Craveling  in  Che  pipe.  Second,  when 
mulciple  cavicies  were  predicCed,  Chey  did  noC  always  physically  exisc. 
Thus,  as  chese  nonphysical  cavicies  opened  and  closed,  addicional  fluc- 
CuaCions  were  inCroduced  inCo  Che  soluCion.  And  Chird,  mathematical 
instabilicies  could  occur  when  Che  gas  bubbles  underwenC  very  large 
volume  changes  during  Che  course  of  a  single  Cime  sCep.  So,  some  cauCion 
is  required  for  Che  choice  of  time  sCeps. 

In  an  efforc  Co  eliminaCe  Che  nonphysical  oscillacions ,  several 
possible  solucions  were  aCCempCed.  The  firsC  "fix"  was  drawn  from  a 
paper  by  Simpson  and  Wylie  (19).  They  showed  thaC  Che  inCegraCion  mechod 
used  when  compuCing  Che  caviCy  volumes  could  inCroduce  some  of  Che  error 
associaCed  wich  nonphysical  cavicies.  In  Che  sCudy  Simpson  and  Wylie 
demonsCraCed  Chree  differenC  meChods  by  which  Eq.  (3.14)  could  be 
integraCed.  The  firsC  was  Euler's  one-sCep  mechod.  This  gave  AVg* 
AC(Q-QU) ,  where  Q  and  Qu  were  Che  flow  upsCream  and  downsCream  of  Che 
caviCy  during  Che  previous  Cime  sCep.  The  nexC  mechod  was  Che  improved 
Euler's  mechod,  which  was  used  in  Chis  sCudy  Co  obCain  Eq.  (3.13).  Here, 
Che  difference  in  Che  flows  aC  Che  previous  and  currenC  cime  sCeps  was 
averaged  Co  find  Che  change  in  volume.  The  Chird  and  lasC  was  Che  for¬ 
ward  inCegraCion  mechod,  which  gave  AVg*AC (Qp-Qpu) ,  where  Qp  and  Qpu  were 
Che  flows  aC  Che  currenC  Cime  sCep. 


The  results  produced  by  using  the  improved  Euler's  method  have 
already  been  presented  in  Figs.  5-8  and  5-9.  Plots  of  the  results 
obtained  by  the  other  two  integration  methods  appear  in  Figs.  5-15  and 
5-16.  A  comparison  of  the  results  in  Figs.  5-8,  5-15,  and  5-16  with 
Fig.  5-7  showed  that  for  this  case  the  improved  Euler's  method  best 
reproduced  the  timing  of  the  main  surge  when  only  vapor  was  accounted 
for.  Both  of  the  other  methods  underestimated  the  cavity  volumes,  so 
that  the  higher  effective  wave  speeds  caused  the  surge  to  occur  sooner 
than  expected.  The  Euler's  method  did  succeed  in  eliminating  the  fluc¬ 
tuations,  but  sacrificed  any  claim  to  accuracy  concerning  the  timing  or 
the  magnitude  of  the  pressure  surge. 

Wylie  (17)  developed  another  method  to  reduce  the  nonphysical 
spikes,  an  air  release  model,  which  was  presented  under  Cavitation  in 
section  III,  Boundary  Conditions.  After  incorporating  Wylie's  air 
release  model  into  the  program  developed  in  this  study,  the  first  simula¬ 
tion  was  repeated.  The  results  of  this  experiment,  at  x=813  m  only,  are 
shown  in  Fig.  5-17.  Upon  comparison  of  Figs.  5-8  and  5-17,  it  was  seen 
that  the  air  release  model  did  reduce  the  magnitude  of  the  oscillations 
to  a  small  degree,  but  did  not  eliminate  them.  Further,  the  air  release 
model  predicted  the  pressure  surge  would  occur  later  than  indicated  in 
Fig.  5-7.  The  small  delay  came  because  the  air  release  model  slightly 
overestimated  the  amount  of  gas  released.  The  additional  gas  was  due  to 
the  inclusion  of  a  very  small  amount  of  air  release  from  the  fluid  along 
with  the  usual  vapor  formation. 

The  last  approach  in  the  attempt  to  eliminate  the  numerical  oscilla¬ 
tions  was  to  reduce  the  size  of  the  time  step.  During  all  the  previous 
simulations  a  time  step  of  0.05  s  was  used.  For  the  results  seen  in  Fig. 


Fig.  5-18.  Results  for  air  release  end  reduced  time  step  (x*813  m, 
L«3048  m,  a*981  m/s,  D“0.61  m,  Qo-0.89  s,  Hv*-10.06  m) 
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5-18,  the  time  step  was  decreased  to  0.01  s.  Unfortunately,  it  was  not 
possible  to  futher  reduce  the  time  step  significantly.  The  number  of 
large  arrays  required  to  complete  the  computation  with  a  smaller  time 
step  overflowed  the  available  memory  in  the  computer.  However,  the  trend 
observed  by  comparing  Figs.  5-17  and  5-18,  indicated  that  the  magnitude 
of  the  oscillations  was  decreased  by  reducing  the  step  size.  It  must  be 
remembered  though,  that  by  reducing  the  step  size,  the  run  time  and  the 
memory  required  were  increased,  so  a  tradeoff  between  accuracy  and  cost 
must  be  made. 


Multiple  Pipe  Verification  of  Valve-in-line 

The  first  of  the  multipie  pipe  simulations  was  taken  from  a  study 
by  Swaf field  (13),  who  was  interested  in  the  flow  downstream  of  a  valve. 

In  essence,  since  the  pipes  of  Swaf field's  system  all  had  the  same 
diameter,  the  region  downstream  of  the  valve  was  analogous  to  the  system 
in  the  first  cavitation  simulation.  The  disturbance  at  the  valve  was 
propagated  downstream,  where  it  was  reflected  in  the  opposite  sense  from 
the  constant  pressure  reservoir  at  the  end.  The  reflection  traveled 
back  upstream  to  the  closed  valve,  where  it  was  reflected  in  the  same 
sense,  propagating  once  more  downstream. 

First  Simulation.  As  seen  in  section  IV,  the  initial  conditions 
for  the  first  of  the  three  cases  were:  Pd“102  kN/m^  abs ,  Qo=*0. 00355  m^/s, 
and  the  valve  cut-off  time  was  0.081  s.  The  results  published  by 
Swaffield  for  this  case  appear  in  Fig.  5-19,  while  the  results  produced 
by  the  program  can  be  seen  in  Fig.  5-20.  The  program  accurately  pre¬ 
dicted  the  magnitude  of  the  first  peak,  but  predicted  magnitudes  which 
were  higher  than  expected  for  the  second  and  third  peaks.  Furthermore, 
the  pulses  predicted  by  the  program,  Fig.  5-20,  occurred  between  the 


pulses  predicted  in  the  two  sets  of  Swaf field's  results  in  Fig.  5-19. 

All  the  results  agreed  that  as  the  valve  closed,  a  cavity  formed  on 
the  downstream  side  of  the  valve.  The  discrepancies  past  the  first  peak 
resulted  from  different  estimates  of  the  amount  of  gas  generated  during 
the  formation  of  the  initial  cavity.  The  disagreement  about  the  cavity 
volume  was  inferred  from  the  differences  in  the  periods  of  the  oscilla¬ 
tions.  The  amount  of  gas  predicted  by  the  program  was  less  than  the 
observed  amount,  but  greater  than  the  amount  predicted  for  pure  vapor. 
Therefore,  the  effective  wave  speed  fell  in  between  the  two  predictions 
by  Swaffield,  causing  the  pulses  also  to  fall  in  between  the  times 
predicted  in  Fig.  5-19.  In  addition,  since  the  cavity  implied  by  the 
program's  results  was  smaller  than  the  cavity  implied  from  the  observed 
results,  the  smaller  cavity  did  not  provide  sufficient  damping  for  the 
later  pulses.  Consequently,  there  was  a  digression  in  magnitude  as  well 
as  period. 

Second  Simulation.  For  the  second  of  the  valve-in-line  simula¬ 
tions,  the  downstream  pressure  was  raised  to  222  kN/m^  abs  and  the  valve 
cut-off  time  was  0.16  s.  A  comparison  of  the  published  results.  Fig. 
5-21,  with  the  program  results  seen  in  Fig.  5-22  and  5-23,  showed 
favorable  agreement  as  far  as  the  general  shapes  and  magnitudes  of  the 
curves  were  concerned.  However,  the  characteristic  closure  curve  used 
in  the  program  simulation  was  not  completely  accurate.  The  approximate 
curve  had  a  more  sudden  closure,  creating  a  slightly  stronger  expansion. 
The  difference  in  the  closure  curves  was  seen  in  the  rate  at  which  the 
pressure  initially  dropped  to  vapor  pressure.  As  a  result  of  this  clo¬ 
sure  the  program  predicted  the  peak  pressure  to  occur  at  0.24  s,  about 
0.1  s  after  the  observed  peak  occurred.  This  delay  implied  a  slight 


overestimate  of  the  gas  released,  causing  a  lover  effective  vave  speed. 
Also,  it  was  possible  to  determine  that  the  gas  release  was  restricted  to 
the  region  immediately  downstream  of  the  valve.  This  observation  stemmed 
from  a  comparison  of  the  results  at  4.04  m  downstream  of  the  valve  with 
the  results  at  S0.8  mm  downstream  of  the  valve.  First,  the  pressure 
4.04  m  downstream  of  the  valve  never  reached  vapor  pressure  during  the 
time  preceeding  the  first  pressure  peak.  And  second,  the  pressure 
oscillations  between  the  first  and  second  large  peak  in  the  results  at 
4.04  m  had  a  period  of  0.02  s,  very  nearly  2L/a  for  the  combined  length 
of  the  second  and  third  pipes  of  the  system,  even  though  vapor  pressure 
was  indicated  at  the  valve  during  this  same  time  interval. 
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Fig.  5-21.  Swaffield's  results  for  case  2  (1*2, 3*9.87  m,  a»918  m/s, 
D»0.0508  m,  Pv»700  Pa,  tc*0.16  s) 

Third  Simulation.  For  the  third  case,  the  initial  conditions  were: 
P^«120  kN/m^  mbs ,  Qo*0. 00541  m^/s,  and  the  valve  cut-off  time  was  0.14  s. 
Once  again,  the  agreement  between  the  published  results  in  Fig.  5-24  and 
the  program  results  in  Fig.  5-25  was  good.  At  50.8  mm  downstream  of 
the  valve,  the  program  predicted  a  peak  pressure  of  approximately  1600  kPa 
at  0.49  s,  which  was  very  close  to  the  results  observed  by  Swaffield. 
However,  in  this  simulation  only  the  results  up  until  the  first  peak 
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were  presented  by  Swaffield.  As  in  the  previous  simulations  the 
agreement  was  good  up  to  that  first  peak.  Unfortunately,  none  of  the 
subsequent  peaks  were  presented  for  this  case.  So,  it  was  impossible  to 
determine  if  the  program  results  for  this  case  were  any  better  than  the 
previous  cases  for  agreement  after  the  first  peak. 

Verification  of  Multiple  Pipe  Valve  Closure 

The  second  of  the  multiple  pipe  simulations  was  also  the  second  of 
the  two  valve  closure  problems  taken  from  Wylie  and  Streeter  (14:60). 

Like  the  previous  multipipe  simulations,  this  system  was  composed  of 
three  pipes  in  series.  For  convenience  the  pipes  were  labelled  as  pipes 
1,  2,  and  3  starting  at  the  upstream  end  of  the  system.  Unlike  the  other 
multipipe  system,  this  one  by  Wylie  and  Streeter  also  contained  area 
changes  at  each  of  the  internal  junctions.  Their  results  for  the 
pressure  at  the  valve  were  tabulated  along  with  the  results  obtained  from 
the  program  and  presented  in  Table  VIII.  A  point  by  point  comparison  of 
the  two  sets  of  results  showed  a  variation  of  less  than  0.1  percent  for 
the  duration  of  the  simulation.  However,  looking  down  the  columns  of 
Table  VIII  or  at  the  appropriate  curve  in  Fig.  5-26,  one  does  not  observe 
the  smooth  rise  and  fall  in  the  pressure  seen  in  the  single  pipe  results, 
Fig.  5-1,  despite  a  basic  similarity  between  the  two  systems.  As 
expected,  in  Fig.  5-26,  the  area  changes  and  the  different  closure  curve 
superimposed  irregularities  over  results  approaching  those  in  Fig.  5-1. 

For  the  first  0.6  s,  the  valve  closed  rather  rapidly  as  the  valve 
characteristic,  r,  dropped  from  1.0  to  0.2.  In  turn  the  hard  closure 
during  this  time  generated  increasingly  stronger  compression  waves,  since 
the  greater  diameter  columns  of  water  in  pipes  1  and  2  had  greater  inertia 
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Table  VIII.  Head  at  the  Valve  for  Multiple  Pipe  Closure 


Time  (s) 

Wylie’  s  Calc. 
Head  (m) 

Program  Calc 
Head  (m) 

0.0 

100.00 

99.99 

0.1 

1-27.65 

127.64 

0.2 

167.51 

167.50 

0.3 

224.67 

224.64 

0.4 

311.71 

311.64 

0.5 

448.71 

448.57 

0.6 

668.70 

668.39 

0.7 

673.58 

673.24 

0.8 

651.84 

651.50 

0.9 

690.25 

689.89 

1.0 

736.11 

735.72 

1.1 

764.86 

764.43 

1.2 

790.15 

789.68 

1.3 

805.23 

804.72 

1.4 

805.76 

805.21 

1.5 

773.  19 

772.64 

1.6 

684.02 

683.51 

1.7 

683.85 

683.32 

1.8 

686.38 

685.81 

1.9 

570.22 

569.71 

2.0 


407.59 


407.18 


Chan  the  column  in  pipe  3.  Aside  from  affecting  the  compression,  the 
differences  in  cross  sectional  area  had  a  second  effect.  At  each  junc¬ 
tion  where  a  pipe  was  joined  to  a  second  pipe  of  different  diameter,  a 
traveling  wave  was  both  transmitted  with  the  same  sign  and  reflected 
with  an  opposite  sign.  This  created  multiple  wave  fronts  within  the 
system.  Thus,  at  junction  two  between  pipes  2  and  3,  the  compression 
waves  traveling  upstream  from  the  valve  were  transmitted  without 
changing  sign  and  reflected  as  expansion  waves.  This  first  reflection, 
referred  to  as  the  first  wave  front,  arrived  back  at  the  valve  in  about 
0.2  s,  but  was  only  strong  enough  to  reduce  the  rate  of  the  pressure 
rise.  The  compression  waves  generated  at  the  valve  were  stronger  than 
the  reflected  expansion  waves  due  to  attenuation  from  friction  and  to 
the  relative  difference  in  time  at  which  each  wave  was  generated. 

Starting  at  0.4  s,  compression  waves  due  to  the  reflection  of  the 
first  wave  front  from  junction  two  arrived  at  the  valve.  Accordingly, 
the  pressure  began  to  rise  more  quickly.  At  0.6  s  the  slope  of  the 
valve  characteristic  decreased  such  that  it  took  the  next  1.2  s  for  the 
valve  to  completely  close.  This  softer  closure  generated  weaker 
compression  waves  at  the  valve.  Thus,  the  combination  of  the  expansion 
waves  from  the  first  wave  front  with  those  from  the  second  reflection  of 
the  first  wave  front  was  stronger  than  the  combination  of  the  current 
compression  waves  with  those  of  the  first  reflection  of  the  first  wave 
front.  So,  the  pressure  began  to  drop  as  the  flow  began  to  reverse. 
However,  around  0.8  s,  the  third  reflection  of  the  first  wave  front 
arrived  about  the  time  the  expansion  waves  resulting  from  the  hard  clo¬ 
sure  in  the  first  wave  front  came  to  an  end.  Therefore,  the  compression 
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Fig.  5-26.  Program  results  for  case  1  of  multipipe  closure  (Li*351  m, 
L2“483  m,  L3-115  m,  a-1200  m/s,  H0-289  m) 

waves  dominated  once  more  and  the  pressure  began  to  rise  again. 

Beginning  at  1.0  s  a  new  factor  entered  the  picture  as  a  second 
wave  front,  consisting  of  expansion  waves  due  to  the  reflection  of  the 
original  compression  waves  from  the  junction  1  between  pipes  1  and  2, 
arrived  at  the  valve.  The  arrival  of  these  waves  initially  caused  the 
rate  of  the  pressure  to  slow.  As  the  strength  of  the  expansion  waves 
increased,  the  rate  of  the  pressure  rise  decreased  even  further.  But, 
before  the  pressure  could  start  to  drop,  compression  waves  due  to  the 
first  reflection  from  junction  2  of  the  second  wave  front  arrived  at  the 
valve  around  the  1.2  s  point.  These  new  waves  acted  to  partially  negate 
the  expansion  waves,  slowing  the  reversal  of  the  flow.  However,  at  1.4 
s,  expansion  waves  produced  by  the  second  reflection  of  the  second  wave 
front  from  junction  two  appeared  at  the  valve  causing  the  pressure  to 
finally  begin  decreasing. 
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It  was  at  1.6  s  that  the  third  wave  front,  created  by  the  original 
compression  waves  reflecting  from  the  reservoir,  reached  the  valve.  At 
this  same  time  two  other  events  occurred.  First,  the  expansion  waves  in 
the  first  reflection  of  the  second  wave  front  due  to  the  hard  closure 
came  to  an  end.  And  second,  compression  waves  due  to  the  third  reflec¬ 
tion  of  the  second  wave  front  arrived  at  the  valve.  These  three  events 
interacted  with  the  other  waves  at  the  valve  to  produce  the  nearly 
constant  pressure  seen  between  1.6  and  1.8  s. 

Another  complex  interaction  occurred  at  1.8  s.  First,  at  this  time 
the  valve  was  completely  closed,  so  no  new  compression  waves  were  being 
generated.  Second,  the  compression  waves  from  the  hard  closure  in  the 
first  reflection  of  the  second  wave  front  ended,  leaving  the  weaker 
waves  from  the  soft  closure.  Third,  compression  waves  due  to  the  first 
reflection  of  the  third  wave  front  reached  tie  valve.  And  fourth,  expan¬ 
sion  waves  due  to  the  fourth  reflection  of  the  second  wave  front  arrived 
at  the  valve.  The  end  result  of  all  these  events  was  a  steep  pressure 
drop. 

Fig.  5-26  includes  a  plot  of  the  results  from  Table  VIII  as  well  as 
the  results  obtained  for  the  simulation  adapted  as  described  in  section 
IV  for  RP-1  and  LH2«  For  all  three  fluids  the  timing  of  the  waves  was 
very  similar  since  the  average  of  the  adjusted  wave  speed  for  all  three 
were  nearly  equal.  The  wave  speeds  for  water  and  for  LH2  were  both 
calculated  as  having  an  average  wave  speed  of  1175  m/s ,  while  RP-1  was 
found  to  have  a  wave  speed  of  1125  m/ s .  Therefore,  since  even  the  wave 
speed  stayed  about  the  same  despite  changes  in  the  fluids,  the  differen¬ 
ces  observed  in  the  transients  in  Fig.  5-26  were  due  to  the  only  remaining 
variables  which  could  change,  either  the  density  or  the  friction  factor. 


In  the  single  pipe  closure  problem,  the  peak  pressure  decreased  by 
about  7  percent  when  the  friction  factor  was  doubled.  In  the  results 
for  the  multipipe  closure,  the  peak  pressure  for  RP-1  was  14  percent 
lower  than  that  of  water,  even  though  the  friction  increased  by  a  factor 
of  1.25.  It  would  be  very  easy  at  this  point  to  reason  that  the  change 
in  density  was  the  major  driver  for  the  significant  change  in  peak 
pressure.  Unfortunately,  it  was  not  that  clear  cut.  The  multipipe 
system  was  longer  and  narrower  than  the  single  pipe  system,  increasing 
the  frictional  effects.  So,  it  was  difficult  to  determine  the  true 
contribution  of  either  the  friction  factor  or  the  density. 

When  the  results  using  LH2  were  compared  with  those  for  water,  a 
much  clearer  indication  of  the  effects  of  density  and  friction  factor 
was  attained.  For  the  very  high  friction  factor  case  for  LH2,  the  atte¬ 
nuation  of  the  reflected  expansion  waves  was  very  great.  Thus,  the 
compression  produced  by  the  valve  dominated  until  the  valve  was  almost 
completely  closed,  only  then  did  the  pressure  peak.  The  situation  was 
very  different  for  the  low  friction  factor  case.  The  pressure  peaked 
relatively  quickly  due  to  the  strong  influence  of  the  reflected  waves  on 
the  pressure  at  the  valve.  But  even  though  the  time  of  the  peak 
pressure  was  very  different,  the  magnitudes  of  the  peak  pressures  dif¬ 
fered  by  less  than  10  percent.  This  small  change  was  in  contrast  to  the 
difference  of  over  100  percent  between  the  peak  pressure  experienced  for 
water  and  the  peak  pressure  for  the  low  friction  LH2.  This  significant 
effect  of  density  was  as  expected  when  the  relationship  of  density  and 
inertia  was  considered.  As  the  density  increased,  the  mass  of  fluid 
within  the  pipe  at  any  specified  time  increased,  which  in  turn  increased 
the  inertia  of  the  fluid.  Thus,  as  the  inertia  increased,  the  force 


required  Co  bring  chat  mass  to  rest  increased,  thereby  increasing  the 
force  per  area  which  can  be  defined  as  pressure. 

Effects  of  Reservoir  Pressure.  The  discussion  above  illustrated 
the  influence  of  density  on  the  severity  of  the  transients  by  a  com¬ 
parison  of  the  results  for  the  different  fluids  in  Fig.  5-26.  In  a 
similar  manner,  a  comparison  of  the  results  for  the  same  fluid  in  Figs. 
5-26,  5-27,  and  5-28  provided  insight  into  the  effect  of  different 
reservoir  pressures  on  the  transient  pressures  in  a  multipipe  system. 

The  first  and  most  obvious  effect  of  changing  the  reservoir 
pressure  was  the  corresponding  change  in  the  initial  pressure  at  the 
valve.  Given  that  the  friction  remained  constant  for  any  particuar 
fluid  in  these  simulations,  the  loss  in  pressure  between  the  reservoir 
and  the  valve  was  always  a  constant  for  a  particular  flow  rate.  So,  as 
the  initial  upstream  pressure  changed,  the  initial  valve  pressure  also 
changed. 

A  second  effect  of  changing  the  reservoir  pressure  appeared  to  be  a 
change  in  the  strength  of  the  reflections.  It  was  expected  that  the 
strength  of  the  waves  generated  during  the  different  closures  would  vary 
as  a  function  of  the  reservoir  pressure.  However,  a  comparison  of  the 
results  seen  in  Figs.  5-26,  5-27,  and  5-28  indicated  that  the  strength 
of  the  reflections  was  also  a  function  of  reservoir  pressure.  The 
results  for  the  different  fluids  retained  the  same  relationship  with 
respect  to  the  results  for  water  as  seen  in  the  original  case  for  this 
simulation.  Therefore,  the  discussion  of  the  second  effect  will  be 
restricted  only  to  the  results  obtained  for  water. 

For  the  reservoir  pressure  of  the  original  simulation,  the  results 
contained  essentially  three  peaks,  as  seen  in  Fig.  5-26.  The  first  peak 
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occurred  when  the  closure  changed  from  hard  to  soft.  A  second  peak, 
also  the  highest,  occurred  just  prior  to  the  arrival  of  the  third  wave 
front.  Lastly,  a  third  peak  that  occurred  about  the  time  the  valve 
finished  closing.  In  the  case  of  the  higher  reservoir  pressure,  Fig. 
5-27,  the  first  peak  occurred  much  the  same  as  in  the  original  case. 
However,  the  second  and  third  peaks  were  not  simply  higher  magnitude 
versions  of  the  peaks  in  Fig.  5-26.  The  second  peak  reached  its  maximum 
as  soon  as  the  second  wave  front  arrived,  an  indication  that  the 
strength  of  the  reflected  expansion  waves  was  enhanced.  A  further 
demonstration  of  the  increased  strength  of  the  reflections  was  that  the 
third  peak  did  not  truly  exist.  There  was  a  leveling  of  the  curve 
around  the  point  where  the  third  peak  should  have  been,  but  the  results 
were  so  dominated  by  the  reflected  expansion  waves  that  there  was  only 
this  change  in  the  slope,  not  a  real  peak.  In  the  case  of  the  lower 
reservoir  pressure,  there  was  only  one  true  peak  remaining  in  the 
results  in  Fig.  5-28.  Unlike  the  high  pressure  case,  the  absence  of 
three  distinct  peaks  was  due  to  the  weakness  of  the  reflected  expansion 
waves,  they  were  unable  to  significantly  affect  the  compression  at  the 
valve.  Thus,  when  each  of  the  first  two  wave  fronts  arrived  at  the 
valve,  there  were  only  changes  in  the  slope.  The  only  true  peak  did  not 
occur  until  the  valve  was  completely  closed.  At  this  time  compression 
waves  were  no  longer  generated  at  the  valve,  so  the  expansion  waves  were 
now  sufficient  to  cause  a  drop  in  the  pressure. 

Effects  of  Valve  Characteristic  Curve.  Moving  now  to  the  results 
seen  in  Fig.  5-29  and  5-30,  the  effect  obtained  by  changing  the  shape  of 
the  closure  curve  was  examined  for  the  multipipe  closure.  As  before, 
the  comparison  of  the  results  concentrated  on  the  results  for  water. 


Starting  with  the  harder  closure,  a  comparison  of  Fig.  5-26  and  5-29 
indicated  two  major  differences  between  the  results  for  the  hard  closure 
and  those  for  the  original  closure.  First,  since  the  compression  waves 
were  stronger  during  the  more  sudden  closure,  the  first  two  peaks  seen  in 
the  results  in  Fig.  5-29  were  higher  in  magnitude  than  the  corresponding 
peaks  in  Fig.  5-26.  This  result  paralleled  the  findings  of  Fig.  5-4 
which  was  a  comparison  of  different  cut-off  times  for  the  single  pipe 
closure.  And  second,  since  the  valve  was  completely  closed  0.6  s 
earlier  in  the  hard  closure  case,  a  greater  pressure  drop  occurred  after 
the  second  peak  than  seen  in  the  original  closure.  In  addition,  the 
change  in  slope  at  the  break  point  of  the  closure  curve  was  more  severe 
for  the  hard  closure.  Therefore,  when  the  expansion  waves  of  the  third 
wave  front  arrived  at  1.6  s,  the  relative  strengths  of  the  waves  were 
such  that  a  third  peak  still  occurred,  despite  the  absence  of 
compression  waves  generated  by  the  valve. 

As  for  the  soft  closure  case,  Fig.  5-30,  the  changes  in  the  slope 
of  t  were  more  gradual.  During  the  first  0.6  s,  the  pressure  rose  in  a 
manner  similar  to  that  seen  in  Fig.  5-26.  At  the  0.6  s  point  was  the 
first  change  in  the  slope  of  the  closure  curve.  Since  the  changes  in 
slope  were  more  gradual  in  the  soft  closure,  the  actual  slope  of  the 
closure  curve  was  much  steeper  than  the  corresponding  portion  of  the 
curve  in  the  original  closure.  Hence,  the  portion  of  the  closure  curve 
between  0.6  and  1.2  s  was  steeper  than  the  original  curve,  and  so,  the 
compression  waves  generated  at  the  valve  were  stronger  for  this  period 
of  time.  Thus,  the  pressure  continued  to  rise  at  a  rate  which  was 
reflective  of  the  closure  curve.  Finally,  as  in  the  low  reservoir 
pressure  case,  the  pressure  peaked  only  after  the  valve  completely 


closed,  ending  the  generation  of  new  compression  waves  at  the  valve,  and 
thereby  allowing  the  expansion  waves  to  finally  dominate  the  solution. 


Saturn  V  Feedline 


In  this  final  system  an  in-depth  wave  by  wave  description  of  the 
interactions  was  not  included,  as  it  would  quickly  become  too  complex  to 
follow.  Instead,  a  number  of  runs  were  made  with  only  the  pipes  of  the 
system  in  order  to  establish  a  baseline.  These  runs  included  variations 
on  the  input  parameters  in  much  the  same  manner  as  were  previously 
incorporated  in  the  single  and  multiple  valve  closure  simulations.  The 
results  for  the  F-l  feedline  were  then  analyzed  by  a  comparison  with  the 
trends  established  by  earlier  simulations  with  the  simpler  systems. 

Once  the  baseline  was  determined,  additional  components  were  added  and 
analyzed,  until  finally  all  the  components  were  included  in  the  final 


Baseline  Simulations.  Beginning  the  baseline  portion  of  the  F-l 
feedline  analysis,  the  results  from  the  first  run  on  Table  VII  appear  in 
Fig.  5-31.  The  shape  of  the  curve  was  surprising  until  the  system  was 


looked  at  in  terms  of  the  findings  of  the  earlier  simulations.  First, 
in  comparison  to  the  earlier  simulations,  the  feedline  was  a  low  pressure 
system.  Furthermore,  the  changes  of  diameter  at  the  junctions  were 
small  compared  to  those  in  the  multipipe  closure.  The  combination  of 
these  two  factors  meant  that  any  reflections  from  the  internal  junctions 
would  be  extremely  weak,  and  hence  the  system  would  tend  to  act  as  if  it 
was  a  single  pipe.  And  second,  the  overall  pipe  length  was  such  that 
the  cut-off  time  for  the  valve  was  much  greater  than  the  oscillatory 
period  of  the  system.  So,  the  system  would  be  expected  to  behave  in  a 


manner  reminiscent  of  the  short  pipe  results  seen  in  Fig.  5-6. 

The  expected  results  extrapolated  from  the  findings  of  the  earlier 
simulations  agreed  well  with  the  results  in  Fig.  5-31.  The  oscillations 
once  the  valve  was  closed  were  a  smooth  sinusoid  with  a  period  of  0.182 
s,  about  the  same  as  would  be  observed  in  a  single  pipe  having  the  same 
period  as  the  feedline.  And,  as  in  the  short  pipe  results  for  the 
single  pipe,  the  pressure  at  the  valve  was  dominated  by  the  relative 
strengths  of  the  compression  generated  by  the  valve  and  the  expansion 
due  to  the  reflection  from  the  reservoir,  creating  the  large  first  peak. 
Also,  as  per  the  description  of  the  single  pipe  valve  closure,  the  expo¬ 
nential  law  of  the  closure  curve  was  reflected  in  the  shape  of  the  first 
peak,  specifically,  in  the  slope  of  the  curve  as  the  pressure  decreased 
just  after  passing  the  peak  magnitude.  This  lasted  up  until  the  time 
when  the  valve  was  fully  closed,  at  which  time  the  pressure  became  domi¬ 
nated  by  the  reflections,  significantly  dropping  the  pressure. 

For  the  results  of  the  second  run  which  appear  in  Fig.  5-32,  the 
only  change  from  the  first  run  was  that  the  cut-off  time  had  been 
increased  to  3.0  s.  As  was  expected  due  to  Fig.  5-4,  the  magnitudes  of 
both  the  primary  and  the  secondary  peaks  decreased  as  a  result  of  the 
softer  closure  produced  by  an  extended  cut-off  time.  Once  the  valve  was 
closed  the  period  of  the  pressure  oscillations  was  still  0.182  s.  And 
finally,  the  changes  in  the  slope  of  the  closure  curve  were  more  gradual 
for  the  longer  cut-off  time.  Consequently,  the  time  was  longer  until 
the  reflected  expansion  waves  were  strong  enough,  relative  to  the  com¬ 
pression  waves  at  the  valve,  to  begin  decelerating  the  flow  and  thereby 
dropping  the  pressure.  Thus,  the  pressure  rose  gradually  to  a  peak  and 
then  gradually  decreased,  as  the  relative  strengths  of  the  waves  slowly 
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reversed  from  the  compression  being  stronger  to  the  expansion  being 
stronger  during  the  first  peak. 

In  the  third  run,  Fig.  5-33,  the  cut-off  time  was  once  more  1.5  s, 
but  E^j  from  Eq.  (4.1)  was  changed  to  1.0.  This  implied  that  the  closure 
curve  was  now  linear.  Therefore,  once  the  initial  reflection  returned 
from  the  reservoir,  the  rate  of  increase  in  pressure  was  reduced  since 
the  compression  waves  and  their  reflections  were  nearly  equal  in  magni¬ 
tude.  However,  because  the  closure  curve  was  linear,  the  closure  did 
not  soften  as  the  cut-off  time  approached  as  in  the  first  two  runs. 
Instead,  the  flow  ended  abruptly  as  indicated  by  the  sharp  drop  in 
pressure  seen  in  Fig.  5-33  at  tc=1.5  s.  Further,  this  sharp  cut-off 
also  affected  the  pressure  oscillations  after  the  valve  closure.  Even 
though  the  oscillations  had  the  same  period,  the  peak  magnitudes  were 
nearly  equal  to  the  magnitude  of  the  initial  peak.  And  lastly,  when  £„, 
was  1.2,  Fig.  5-31,  the  slope  of  the  closure  curve  started  off  steeper 
than  the  linear  closure,  then  gradually  decreased  until  the  slope  was 
less  than  that  for  Em-1-  Thus,  when  Fig.  5-31  and  5-33  are  compared, 
the  pressure  in  Fig,  5-31  increased  quicker  due  to  the  steeper  slope  of 
the  closure  curve.  As  the  first  reflection  arrived  the  pressure  rise 
slowed  in  each  case.  However,  in  Fig.  5-31,  the  strength  of  the 
compression  decreased  as  the  slope  of  the  closure  decreased,  causing  the 
peak  and  gradual  decline  in  pressure.  While  in  Fig.  5-33,  the  slope  did 
not  change,  so  the  strength  of  tha  compression  remained  the  same. 
Therefore,  since  the  reflections  had  been  attenuated  by  friction,  they 
were  not  as  strong  as  the  compression  and  the  pressure  continued  to 
rise,  albeit  slowly.  This  rise  continued  until  tca,1.5  s,  when  the  valve 
closed,  ending  the  generation  of  compression  waves  at  the  valve.  Due  to 


this  continual  pressure  rise  in  the  results  in  Fig.  5-33,  the  peak 
pressure  was  slightly  higher  than  that  seen  in  the  results  in  Fig.  5-31. 

The  results  which  appeared  in  Fig.  5-34  were  those  for  run  4,  where 
Em“1.5.  In  this  case  the  closure  was  initially  harder,  but  ended  softer 
than  the  closure  experienced  in  Fig.  5-31,  producing  several  effects. 

First,  the  slope  of  the  initial  closure  was  steeper  than  the  previous 

runs,  hence  the  peak  pressure  was  higher.  Second,  in  Fig.  5-34,  due  to 
the  more  radical  difference  in  the  slopes  of  the  closure  at  the 
beginning  and  at  the  end,  the  relative  strengths  of  the  compression  and 
its  reflections  caused  the  pressure  to  decrease  more  rapidly  than  in 
Fig.  5-31.  And  third,  during  the  final  moments  of  the  closure,  the 
slope  of  the  closure  curve  more  closely  approached  zero  when  Egj-1.5  than 
in  the  earlier  cases.  For  this  reason  the  magnitude  of  the  pressure 
oscillations  was  small,  though  the  period  still  remained  the  same.  This 

was  expected  due  to  the  results  seen  in  Fig.  5-3  and  5-6  for  the  single 

pipe.  From  Fig.  5-3  it  was  determined  that  the  pressure  rose  more 
quickly  and  peaked  higher  as  Eg,  increased.  And,  except  for  the  fact 
that  the  cut-off  time  was  decreased  to  1.5  s,  Fig.  5-34  could  almost  be 
one  of  the  curves  in  Fig.  5-6,  the  results  agree  that  well. 

Run  5  was  the  last  of  the  baseline  cases.  For  this  case  the  inputs 
had  been  returned  to  the  original  values,  except  for  the  initial  valve 
pressure,  which  was  increased  by  40  psi,  2.76x10^  Pa  in  metric  units. 

The  results,  which  appear  in  Fig.  5-35,  were  almost  an  exact  duplicate 
of  those  in  Fig.  5-31.  The  only  observable  differences  were  the 
increases  in  magnitude  of  the  primary  and  secondary  peak  pressures.  The 
shape  of  the  first  peak  and  the  frequency  of  the  pressure  oscillations 
in  Fig.  5-35  remained  the  same  as  in  Fig.  5-31. 


Fig.  5-35.  Results  of  run  5  for  Saturn  V  feedline  (PQ*896  kPa,  tc-1.5  s 


Results  of  run  6  for  Saturn  V  feedline  (P0»621  kPa,  tc*1.5 
Em«1.2,  w-2.5  Hz) 
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Pulser  Valve  Simulations.  The  next  aix  runs  began  the  inclusion  of 
other  components  into  the  simulations.  For  these  runs  a  pulser  valve  was 
included,  operating  at  frequencies  of  2.5,  10,  and  25  Hz.  The  first  and 
third  frequencies  were  chosen  to  test  the  frequency  response  of  the 
models  in  use,  the  frequencies  were  arbitrarily  chosen  near  two  resonant 
frequencies  reported  for  this  system  by  Dehoff  (2).  The  last  frequency 
was  chosen  as  a  control  to  provide  a  median  value  which  was  not  iden¬ 
tified  as  a  resonant  frequency. 

The  results  for  the  first  of  these  runs  appeared  in  Fig.  5-36.  In 
this  run  the  frequency  of  the  pulser  was  2.5  Hz.  The  results  seem  to 
indicate  a  resonant  condition  due  to  the  large  oscillations  that  built 
up  and  finally  masked  the  expected  response.  It  should  be  noted,  that 
due  to  the  choice  of  the  function  to  represent  the  flow  out  of  the 
pulser  valve,  Eq.  (4.5),  the  frequency  of  the  sinusoid  was  effectively 
doubled  due  to  the  absolute  value  being  used.  Hence,  the  frequency  of 
2.5  Hz  became  5  Hz  for  all  practical  purposes,  and  very  nearly  the 
nearly  the  natural  frequency  observed  in  the  previous  runs.  Thus,  the 
system  was  near  resonance,  ana  the  magnitude  of  the  pressure  oscilla¬ 
tions  became  so  large  that  the  system  began  to  reach  vapor  pressure. 
Consequently,  as  with  the  cavitation  simulations,  pressure  spikes  were 
introduced  into  the  solution.  More  important  though,  the  oscillations 
remained  bounded  and  finally  decreased  because  the  gas  released  by  the 
column  seperation  changed  the  effective  wave  speed,  thereby  changing  the 
natural  frequency  of  the  pipeline. 

When  the  frequency  of  the  pulser  was  increased  to  10  Hz  as  in  Fig. 
5-37,  the  results  were  not  as  dramatic  as  those  seen  in  Fig.  5-36.  In 
this  case  the  results  resembled  Fig.  5-31  with  a  small  magnitude 
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37.  Results  of  run  7  for  Saturn  V  feedline  (Pa-621  kPa, 
Ea«1.2,  «-10  Hz) 


oscillation  superimposed  on  top  of  them.  At  this  frequency  the  oscilla¬ 
tions  at  the  valve  started  small  and  built  up  to  a  magnitude  that  was 
only  a  fraction  of  the  magnitude  of  the  oscillations  following  the  valve 
closure.  Increasing  the  frequency  further  to  25  Hz  produced  the  results 
seen  in  Fig.  5-38.  Once  again,  the  oscillations  were  simply  superim¬ 
posed  over  the  results  of  Fig.  5-31.  The  pulser's  effective  frequency 
of  50  Hz  was  nearly  a  harmonic.  This  was  indicated  in  the  results  by 
the  magnitude  of  the  oscillations  growing  more  quickly  and  to  a  larger 
value  than  those  seen  in  Fig.  5-37. 

The  last  three  runs  of  this  set  examined  the  effect  of  changes  in 
some  of  the  system  parameters  while  the  pulser  was  present.  In  run  9, 

Eg  was  changed  to  1.0  to  see  if  a  change  in  the  closure  curve  affected 
tne  resonance  condition  exhibited  in  Fig.  5-36.  As  seen  by  the  results 
in  Fig.  5-39,  the  change  had  no  significant  impact.  Jumping  from  the 
case  with  the  largest  oscillations  to  the  one  with  the  smallest,  for  run 
10,  run  7  was  repeated  with  the  initial  valve  pressure  increased  to  130 
psi  which  is  8.96x10^  Pa.  Examining  the  results  of  run  10  in  Fig.  5-40, 
the  curve  was  little  more  than  a  translation  in  the  y-direction  of  Fig. 
5-37.  However,  as  seen  by  Fig.  5-41,  when  the  system  pressure  was 
increased  while  keeping  Eg-1.0,  the  resonant  condition  was  not  as  domi¬ 
nant.  While  still  very  sizable,  in  this  case  the  pressure  oscillations 
after  the  valve  closed  were  out  of  phase  with  the  oscillations  due  to 
the  pulser.  Hence,  the  erratic  pressure  response  after  the  valve  was 
fully  closed. 

Accumulators  and  Cavitation  Simulations.  The  next  six  runs  incor¬ 
porated  the  other  system  components,  the  "pump  inlet  cavitation"  and  the 
prevalve  accumulator.  Run  12  started  things  off  with  cavitation  only. 


Upon  comparison  of  the  run  12  results,  Fig.  5-42,  with  the  results  in 
Fig.  5-33,  one  major  change  could  be  noted.  Once  the  valve  was  closed, 
the  results  in  Fig.  5-42  showed  a  decrease  in  the  frequency  from  5.5  Hz 
to  3.2  Hz.  This  decrease  in  frequency  was  to  be  expected  however.  In 
run  12  there  was  now  vapor  present  in  the  form  of  the  turbopump  inlet 
cavitation  which  produced  a  reduction  of  the  effective  wave  speed  in 
that  region,  thereby  increasing  the  travel  time  of  the  waves. 

In  runs  13,  14,  and  15  the  cavitation  was  still  present,  however, 
the  prevalve  accumulator  had  been  included  with  volumes  as  seen  in  Table 
VII.  Comparing  Fig.  5-42  with  Figs.  5-43,  5-44,  and  5-45  deliniated 
several  trends.  The  most  notable  trend  was  the  increase  in  the  period 
of  the  pressure  oscillations  as  the  gas  volume  in  the  accumulator 
increased.  Also,  the  strength  of  the  reflections  from  an  internal  boun¬ 
dary  condition  increased  with  increasing  gas  volume.  This  aspect  was 
best  demonstrated  by  the  increase  in  the  rippling  of  the  pressure  trace 
as  the  accumulator  gas  volume  increased.  And  finally,  there  was  no 
significant  decrease  in  the  peak  magnitudes  as  the  accumulator  volume 
increased.  Taken  together,  these  trends  indicated  a  problem  in  the 
modeling  of  an  accumulator.  Instead  of  acting  as  a  surge  suppressor, 
the  model  developed  for  this  study  acted  in  a  manner  very  similar  to  a 
vapor  cavity  in  the  discrete  bubble  model  of  cavitation. 

This  supposition  concerning  the  accumulator  model  was  further 
demonstrated  by  the  results  of  runs  16  and  17.  In  run  16  the  initial 
valve  pressure  was  again  increased  to  896  kPa .  In  comparison  to  the 
results  in  Fig.  5-42  the  results  were  as  expected.  The  pressure  trace 
was  essentially  the  same  curve  translated  in  the  y-direction.  The  fre¬ 
quency  of  the  pressure  oscillations  even  stayed  the  same.  However,  com- 


Fig.  5-43.  Results  of  run  13  for  Saturn  V  feedline  (P0“621  kPa,  tc“1.5 
Em“1.0,  vo 1. *0.018  m^) 


Fig.  5-44.  Results  of  run  14  for  Saturn  V  feedline  (P0“621  kPa,  tc*1.5 
Ea“ 1.0,  vol. “0.036  «3) 


paring  between  the  results  of  run  16  in  Fig.  5-46  and  those  of  run  17  in 
Fig.  5-47,  the  period  of  the  oscillation  increased  significantly  when 
the  accumulator  with  an  initial  volume  of  0.054  m^  was  included. 
Furthermore,  the  existance  of  relatively  strong  internal  reflections  was 
once  again  readily  apparent.  And,  there  was  still  no  surge  suppression 
demonstrated  by  the  accumulator  model.  Thus,  when  all  the  components 
were  included  for  run  18,  the  massive  pressure  oscillations  seen  in 
Figs.  5-36  and  5-39  were  not  eliminated  from  the  results  seen  in  FLg. 
5-48  by  surge  suppression.  Rather,  the  inclusion  of  the  accumulator 
suppressed  the  surges  by  changing  the  natural  frequency  of  the  system 
such  that  pressure  oscillations  caused  by  the  pulser  took  longer  to 
build  up  to  the  larger  magnitudes. 
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VI.  Conclusions 


la  this  study  s  simple  but  flexible  computer  program  was  developed 
and  used  to  analyze  fluid  transients  in  pipelines.  The  derivation  of 
the  equations  of  motion  used  in  the  numerical  analysis  was  based  on  the 
method  of  characteristics.  Boundary  conditions  were  obtained  from 
either  published  sources  or  developed  mathematically  to  allow  for  con¬ 
ditions  which  commonly  occur  in  propellant  feedline  systems. 

The  results  of  the  program  were  quite  good  when  compared  against 
published  data.  Further,  a  number  of  trends  were  observed  in  the  results 
1)  The  integration  method,  the  inclusion  of  air  release,  and  the  size  of 
the  time  step  all  affected  the  numerical  stability  of  the  discrete 
bubble  method.  2)  The  valve  characteristic,  r,  the  cut-off  time,  the 
friction  factor,  and  the  pipe  length  strongly  influenced  the  peak  magni¬ 
tude  of  the  pressure  transient.  3)  The  density  strongly  influenced  the 
peak  magnitude  of  the  transient  pressure  through  inertial  effects.  And 
4)  The  pressure  of  the  system  affected  the  strength  of  the  reflections 
from  the  internal  boundary  conditions. 

It  is  recommended  that  a  new  accumulator  model  be  found  or  created 
to  replace  the  one  currently  included  in  this  study.  Also,  even  though 
the  results  obtained  through  the  use  of  this  program  were  good,  there  is 
still  room  for  additional  improvements  in  the  modeling  and  procedures. 

One  specific  area  that  might  be  improved  is  that  section  of  the  code 
which  provides  the  initial  solution  for  the  numerical  analysis. 
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Computational  Requirements 

The  program  developed  during  this  investigation  was  written  in 
FORTRAN  77  and  implemented  on  the  ASD  CDC-Cyber  mainframe  computer.  The 
routine  which  is  detailed  later  in  this  Appendix  did  not  require  a  great 
deal  of  computational  power.  In  the  case  of  the  single  pipe  models,  run 
times  averaged  up  to  20  CPU's.  The  Saturn  V  model  usually  required 
around  200  CPU's  due  to  the  small  time  step  required  to  achieve  an 
integer  number  of  reaches  in  the  very  short  pipes. 


Program  Explanation 

Fig.  A-l  is  a  general  flow  diagram  of  the  routine  developed  in  this 
study.  The  program  begins  by  dimensioning  all  of  the  arrays  that  are 
found  in  the  program.  Since  FORTRAN  cannot  dimension  arrays  after 
reading  in  any  data,  the  user  must  be  sure  the  arrays  are  large  enough 
at  the  start.  The  next  part  of  the  program  is  the  reading  of  the  input 
data.  The  user  must  pay  close  attention  to  the  sequence  of  Read  state¬ 
ments  in  this  program  (see  program  listing).  The  first  three  Read  sta¬ 
tements  are  used  only  once,  the  remaining  three  Reads  are  repeated 
within  seperate  loops.  Thus,  the  input  deck  must  contain  sufficient 
data  to  satisfy  the  full  range  of  any  one  Read  statement  for  the  length 
of  its  loop  before  trying  to  input  data  for  the  next  Read. 

After  all  the  data  has  been  input,  if  the  wave  speed  has  been  set 
to  zero,  the  program  will  then  calculate  a  value  by  2q.  (4.2).  Once  the 
wave  speed  is  known  each  pipe  must  be  subdivided  into  reaches  for  com¬ 
putational  purposes.  The  routine  the  program  uses  to  make  this  calcula¬ 
tion  is  taken  from  Wylie  and  Streeter  (14:45).  In  each  pipe  of  the 


system  it  is  required  that 


At-L/aH 


in  which  N  is  an  integer.  However,  for  computational  purposes  the  time 


step  must  remain  constant  for  all  the  pipes  in  the  system.  Wylie  and 


Streeter  argued  that  the  wave  speed  does  not  need  to  be  known  exactly, 


thus  may  be  adjusted  slightly  to  arrive  at  an  integer  number  of  reaches 


in  every  pipe.  In  equation  form 


At*L/a(l+c)N 


(A. 2) 


where  c  is  a  constant  equal  to  the  permissible  variation  in  the  wave 


speed,  which  is  always  less  than  0.15. 


With  the  pipes  now  subdivided  for  computation,  the  finite  dif¬ 


ference  routine  needs  an  initial  condition  from  which  to  start.  In  this 


*a»  •  •  •_*  ’j 


■ 


r 


program  the  steady  state  pressure  at  either  the  upstream  or  the 
downstream  end  must  be  known.  If  the  upstream  pressure  is  input,  the 
program  will  Chen  use  the  steady  state  form  of  the  finite  difference 
equations  to  relate  the  upstream  pressure  to  the  pressure  at  any  point 
in  the  flow.  However,  if  the  downstream  pressure  is  input,  the  program 
must  first  work  backwards  through  the  steady  state  solution  to  find  the 
pressure  at  the  upstream  end.  Once  the  pressure  is  known  at  the 
upstream  end,  the  pressure  may  be  calculated  at  any  point  as  before.  In 
the  program  listing  the  loop  to  work  backwards  through  the  solution  is 

DO  22  etc.,  so  the  user  must  be  careful  to  eliminate  this  loop  if  the 

upstream  pressure  is  input  as  an  initial  condition. 

Once  the  program  has  the  initial  condition  from  which  to  start,  the 
transient  loop  begins.  The  first  part  of  this  loop  is  the  calculation 
of  the  pressure  and  the  volume  flow  rate  at  all  the  interior  points  of 
the  pipes.  After  the  values  at  the  interior  points  have  been  calcu¬ 
lated,  a  loop  to  apply  the  boundary  conditions  to  any  internal  junctions 
begins.  First,  the  flag,  M,  is  checked  to  see  if  there  is  more  than  one 

pipe.  As  shown  in  Fig.  A-l,  if  there  is  only  one  pipe  the  program  skips 

down  to  the  exterior  boundary  conditions.  However,  if  there  is  more 
than  one  pipe,  the  various  flags  are  then  checked  to  determine  which  of 
the  boundary  conditions  exist  at  that  particular  junction.  If  no  other 
boundary  condition  is  specified  the  default  boundary  condition  is  a 
series  connection. 

The  last  section  of  the  program  before  the  next  time  increment 
begins  is  the  exterior  boundary  conditions  of  the  system  being  analyzed. 
These  boundary  conditions  are  not  actually  exterior  to  the  pipe  system, 
simply  the  conditions  at  the  upstream  end  of  the  first  pipe  and  the 


downstream  end  of  the  last  one.  The  user  must  make  sure  the  correct  end 
conditions  are  being  used  since  a  number  of  different  end  conditions 
have  been  left  in  the  program  as  Comment  statements. 


Program  Variables 


B,  FF,  R 


CM,  CP 


DPPO ,  PPO 


wave  speed 


pipe  area 


pipeline  constants,  Eq.  (2.20)  and  (2.21) 


B0,B1,B2,B3  pump  constants,  Eq.  (B.8) 


constants  from  Eq.  (2.18)  and  (2.19) 
CM  for  pipe  2  of  a  branch 
constant  from  Eq.  (3.6)  and  (3.8) 


initial  value  of  CV 


pipe  diameter 

pressure  rise  due  to  pump 


pump  constants 


DT  time  step 

E  Young's  modulus 

EM  exponent  from  Eq.  (4.1) 

F  friction  factor 

G  gravitational  acceleration 

I  counter  to  indicate  which  reach  within  a  pipe 

ICAV  indicator  for  the  existence  of  cavitation  bubbles 

IPR  constant  governing  amount  of  output 

J  counter  to  indicate  which  pipe 

JMAX,  IMAX  location  of  maximum  pressure 
J2  counter  for  branches 


I 

|«  1’ 

“w 


I 

I 

I 


K 


counter  used  with  XPR 


1 

1 

1 


8 

I 

« 


§ 

5 


$ 


2 


l 


r 


KL 

L,  LB 
M 
MB 
N 

NA 

NB 

NP 

NPLS 

NS 

NTA 

NV 

OLDCM,  OLDCP 
P 

PAO 

PB 

PMAX 

PO 

POR 

PP 

PS  I 

PV 

Q 

QP 

QPO 

QPU ,  QU 


bulk  modulus 

flags  for  branches  (see  top  of  program  listing) 
number  of  pipes  in  main  line 
M  ♦  NB 

number  of  reaches  in  a  pipe 
accumulator  flag 
number  of  branches 
pump  flag 
pulser  flag 

last  computation  point  in  a  pipe 

number  of  points  if  TAU  is  not  exponential 

valve  flag 

CM  and  CP  from  previous  time  step 
pressure  from  previous  times  step 
initial  pressure  of  an  accumulator 
barometric  pressure 
maximum  pressure 
initial  pressure 

initial  pressure  at  the  upstream  end 
pressure  for  current  time  increment 
adjustment  factor  for  wave  speed 
vapor  pressure 

volume  flow  rate  at  previous  time  step 
volume  flow  rate  at  current  time  increment 
initial  volume  flow  rate  of  a  pump 

volume  flow  rate  on  upstream  side  of  computing  location 
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I 

I 

I 

I 

I 

1 

QV  volume  flow  race  through  pulser  valve  1 

‘  '  f 

RHO  density 

T  time  | 

TAU  closure  character ist ic ,  Eq.  (3.5) 

« 

TC  cut-off  time  for  valve  at  end  of  the  system  ; 

TCUT  cut-off  time  for  internal  valve 

TH  interpolation  factor  for  nonexponential  TAU 

TMAX  time  limit 

TMX  time  of  maximum  pressure 

TOL  tolerance  for  accumulator  iteration 

TT  pipe  wall  thickness 

TX  exponent  of  Eq.  (4.1)  for  valve  at  end  of  system 

VAO  initial  volume  of  accumulator 

VCAV  volume  of  a  vapor  cavity 

W  frequency  of  pulser  valve 

XL  pipe  length 

interpolation  factor  from  Eq.  (2.24) 


ZETA 


PROGRAM  THESIS 

REAL  A<17>,D<17>,F<17>,XL<17>,B(17>,R<17>,FF<17),ARU7> 

REAL  F(17fl50)»PP<17»150)»Q(17»l50)f0P(17,i50)»TA(30)»ZETA(17) 
REAL  TOUT (17) yEM< 17) »TAU( 17r2) .CV( 17f2>  »CVP< 17»2) »TT( 17) 

REAL  OU( 17(150) »QPU< 17ylS0)yVCAV(17yl50)yE(17)yALPH(17) yKL 
REAL  00U<17,130>,0Q<17,130>,RL<17>,V<17>,CI<17> 

REAL  OLDCP( 17)yOLDCM< 17)yT0L(17) yPAO( 17)yVA0(17) 

INTEGER  NS<17>,N<17>,LBU7),L<l7>,NV<17>,ICAV<17,130>,NA<i7> 


COHMOH/DELTAP/  BOy Bl tB2,»3, DPOyPPOyDPPO 
COMMON/OLD/  OLDCH,OLDCP 


MV  INDICATES  IF  THERE  IS  A  VALVE  AT  INLET  OF  PIPE  J 
NA  INDICATES  IF  THERE  IS  AN  ACCUMULATOR  AT  INLET  OF  J 
L  INDICATES  IF  THERE  IS  A  (RANCH  AT  THE  OUTLET  OF  PIPE  J 
IN  EACH  CASE  NO-O  AND  YES-1 
ND  IS  THE  NUMBER  OF  BRANCHES 
LB  INDICATES  WHICH  JUNCTION  BRANCH  J  OCCURS  AT 
NP  INDICATES  WHICH  JUNCTION  A  PUMP  OCCURS  AT 
NPLS  INDICATES  WHERE  A  PULSING  VALVE  IS  LOCATED 


INPUT  OF  SYSTEM  PARAMETERS 


READ<27,*>Q0,P0yGyDT,PSIyM,THAX,IPR,TCyNByKL,RH0 

READ<27,*)NP,QP0yPP0yDPP0yNTA,B0yBl,B2,B3yPB,PV 

READ (27ft)QV»Wy NPLStTXf PHI »AO»DESIG 

MB-M+NB 

DO  10  I-1,MB 


I 

10 

CONTINUE 

DO  11  J-1,MB 

READ(27y«)EM(J)yTCUT(J)yNV(J)yVA0(J)yT0L(J)yALPH<J) 

11 

CONTINUE 

§ 

DO  12  IMyNTA 

READ(27f t)  TA< I ) 

12 

CONTINUE 

ff 

POR-PO 

PMAX-PO 

QOR-QO 

.*■ 

DO  20  J«1,MB 

X* 

C 

C 

C 

ADJUSTMENT  FOR  INTEGER  NO.  OF  REACHES 

IF  (A(J).EQ.O)  THEN 

A( J)-80RT<  <KL/RHO)/< 1+CKLBD( J) )/<E( J)>TT( J) > ) ) 

ENBIF  y 

NL-XL< J)/< ( 1+PSI )>A< J)tDT) 
N<J)»XL(J>/<  <1-PSI)*A( J)BDT) 
IF  <NL.EQ.N<J))  GOTO  15 
XN-XL(J)/<DT*A<J)) 

N<J)-XN 

IF  <XN-N< J) .GT « >5)  THEN 
N<J)-N<JH1 
ENDIF 


o  o  o 


IF  <N<J).EO.O>  QOTO  14 
A( J)«Xt( J)/(N< J)tDT> 

IS  NS<J)«N<JH1 

UNITE!*, •>'.)■  ',J,'  SECTIONS*  ',NS<J> 

14  IF  (N<  J) .EO.O)  THEN 

UNITE!*,*) 'DT  TOO  SIS  FON  FIFE  ',J 
SOTO  9f 
EMI  F 

CALCULATION  OF  SYSTEM  CONSTANTS 

AN!J)».7S54*D!J>**2 

NL< J)«XL< J)/N! J) 

VU)«NL<  J)*AN!J> 

POST-PO+PD 
Cl!J)»POST*AO*V! J) 

ZETA! J)*N( J)*DT*A! J)/XL! J) 

B<  J)*AN< J)/(RHO*A< J) ) 

R(J)-0*AA ( J ) «DT*S!N!  ALPW ( J) > 

FF( J)«F! J)*DT/!D! J)*AN( J)*2) 

NN-NS(J) 

DO  20  I-1'NN 
VCAW( Jr  I )«0. 

ICAU! J,I)*0 
20  CONTINUE 

DO  22  J*1 »H 

P0N-P0N+FF!J>*N!J)*0P0**2/»U> 

22  CONTINUE 

UNITE!4,*)6,P0N,P0,Q0,DT,THAX,IPN 

INITIAL  CONDITIONS 

P !  1 » 1 ) *PON 
FAINT  t,P! 1,1 ) 

DO  2S  J«1,HB 
NN*NS< J) 

IF  (J.6T.1)  THEN 
P!J,l>«P!J-l,N8!J-i>> 

IF  (NU(J).EQ.l)  THEN 
P!J»1)***P*P< J»1 ) 

ENDIF 

ENOIF 

IF  (NA(J).EQ.l)  THEN 
PA0!J)*P!J,1> 

EMXF 

DO  2S  I*1»NN 
Q!J»I)*QPO 
QU! J,I)*QPO 
00! J,I)*0P0 
OQU!  J»I)»QPO 

IF  !J.EO.NP.AND.I.£Q.i)  THEN 
CALL  DELP! J,A,Q,P,NS,DP> 
P!J,l)-P!J-l,NS!J-imDP 
ENDIF 

P!J,I)»P<J,1)-!I-1)*FF!J>*Q!J,I>**2/D!J> 


Ii*  4.*  «-»  i 


>,»  L.*  i/  Ul 


TTTT 


CONTINUE 

H1-N41 

IF  (NI.NC.N)  THEN 
00  30  J*Nl ,Ni 

P<J,l)-P<Li<J),NS<LIU>) ) 
WNNS(J) 

90  30  I»1,NN 
Q<  J, I)-0. 

GU<J,I)-0. 

P<J,I)-P <J,1> 

CONTINUE 

EN9IF 

TAU<N,2)-TA<1> 

CW<H,2>-QP0»*2/(2«P<M,NS(N)>> 

T«0 

K-0 


V 

o 

00  35  J-1,NI 

?• 

CVP(J,l)-QPO*t2/<0.01*P<J,l>> 

UNITE (4,*)J,A(J) ,F( J) , ZETA(J) 

S 

35 

CONTINUE 

UNITE (6,9) 'DESIG-  ',DE3I0 

40 

UNITE <0,i)T,P(H,N8(H) ),P(14,1),P(8,1),Q(1,1> 

C 

UNITE(6,9) 'HAX  AT  ' ,THX, JMAX,IHAX,PHAX 

i 

DO  45  J»1,NS 

c 

UNITE(6,9) ' J-  ' ,J,P<J,1),Q<J,1),P<J,NS<J> 

c 

UNITE(0,9)Q( J,1 ) ,0( J,NS( J>  > 

$ 

c 

UNITE(6,9) ' J"  ',J,(P(J,I)»I-1,NS(J)) 

v' 

c 

MNITE(0,9) <Q( J, I ) , I-l ,NS<  J) ) 

c 

UNITE<0,9> 

1 

45 

CONTINUE 

9 

C 

c 

r 

TNAN8IENT  LOOP 

\ 

w 

50 

T-T4DT 

IF  (T.GT.THAX)  GOTO  ?? 

K-K+l 

9 

C 

c 

I* 

INTEKION  POINTS 

■  ■ 

L 

DO  01  >1,N9 

V 

IF  (T.GE.TCUTU))  GOTO  51 

TAU<J,l>«<l-T/TCUT<J)>90EH<J> 

GOTO  52 

c: 

>. 

51 

TAU<J»l>-0, 

52 

CV<J,l>«CVP<J,l)9TAU<J, 1)992 

IF  (NN.EQ.l)  GOTO  01 
DO  00  1-2, NN 

CALL  FCN( I , J,Q,0U,F,B,ft,FF,2ETA,CH) 
CALL  FCP<I» J»QfQU»P,i»N»FF»ZETA»CP) 
IF  <ICAg<J,I).EQ.l>GOTQ  53 
QP<J,I)*<CPKN>/2 
PPU,I>«<CP-GP<J,I>>/KJ> 

IF  <PP<J,mP9.LT.PV>  GOTO  53 
GOTO  SO 


,N 

•  V 

.  -  -  *  »N 


.«  it  1» 


B 


I 


s 


33 

C 


57 


C 

C 


CONTINUE 

MNITE<t,t> 'CAVITATION  IN  SECTION  OF  PIPE  ',J 

BI*(PHIi<CH-CP)*VCAV< J» I )/<2*DT)i ( l-PHI )*(OQ< Jr I ) 

-OQU< J»I)> )/<PHII8< J) > 

C4*C1< J)/(2*DT*PMI«B( J) ) 

PPU»I>*<-BH-2*<PV-PB>*SQ*T<(8I62»<PV-PB>>*t2f8*C4>>/4 

PP<J,I)—PB4PV 

ICAV<JfI)*l 

QPU<J,I)»CP-PPU,I)*B<J> 

QP< J»I)*CH+PP< J»I)*B( J) 

VCAV< J» I)*VCAV< J»I )+DTt<l*HIt(QP( J»I)-QPU( J» I) ) 

+■(  1-PHI  )f  (G<  J*  I  )-QU(  J»T)  > ) 

VCAV ( J  *  I ) ■  VC  AV (JrI)+2tDT*(PHI*(QP(J,I) -QPU <  J » I ) ) 
><l-PHI)t(OQ( J»I>-OQU< J»  D) ) 


w 

IF  (VCAV(Jyl).GT.O.)  GOTO  60 
ICAV<J,I)-0 

ft 

VCAV( Jyl)*0. 

5? 

QP(J,I)*(CP6CH)/2 

PP< J»I)*<CP-QP< J»I) >/B( J) 

QPU(J,I)HJP(J,I) 

60 

CONTINUE 

61 

CONTINUE 

s 

C 

e 

C 

PIPE  JUNCTIONS 

C 

IF  (N.EQ.l)  GOTO  73 
J2-H 

DO  70  J-2,N 

IF  <L<J-1).NE.0>  THEN 
J2-J2+1 
ELSE 
J2-J2 
ENDIF 
1*1 

CALL  FCH(IfJfQrQUrPfB,ll»FFrZETA»CN) 
CN1-CH 

CALL  FCN<If J2»0fQUrPrB»RrFFf ZETArCN> 

Jl-J-1 

I-NS(Jl) 

CALL  FCP<X'JltQ'QU»PfBfR'FF,ZETA,CP> 
IF  (T.EQ.DT)  THEN 
OLMFUD-CP 
OLDCM<J)-aU 
ENBIF 

IF  <NV<J)»EQ*1)  GOTO  43 
IF  <NA<J).E0.1>  GOTO  64 
IF  (L(J-l).EQ.l)  GOTO  63 
IF  O.EQ.NP)  GOTO  67 
IF  (J.EQ.NFLS)  GOTO  4? 

PP<  J,1  >»<CF-CH1  >  /  <  B<  JHB  <  J-l  >  > 

QP( J»1)*CH1+B< J)tPP< 1) 

QPU< J»1)*QP< J»l) 

QP< J-1»I)*0P< J»l) 

QFU<J-l,I)*OF(J,l) 

PP(J-l,I)*PP(J,l) 


96 


43 

C 


I 

8 

1 

i 

I 

S 


I 

s 

8 


g 

I 

1 


64 

65 
C 

67 

C 

67 

70 

C 

C 

C 

75 

C 

C 

C 

C 

C 

C. 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

1 

c 

c 

c 


77 

80 


GOTO  70 

CALL  VALVE<J»BfNS,CP»CHl,CV,QP,PP,QPU) 

URITE<tf!) 'VALVE  CALLED' 

GOTO  70 

CALL  ACCUN<J,B,P,0,CP,CM,DT,NS»PA0,VA0,T0L,PP,QP,QPU) 
GOTO  70 

CALL  DftANCH<J,J2,D,NS,CP,CH,CHl,QP,PP,QPU) 
NftITE<»»*)'DftANCH  CALLED' 

GOTO  70 

CALL  DELP< J»A»0fP»M8f DP) 

CALL  PUAP( JfB»CH»CP»DP»M8»0P»PPf0PU) 

HftITE(#,»)'PUHP  CALLED' 

GOTO  70 

CALL  PULSE<J,B,T,M,N8,QV,CP,CH1,PP»QP,QPU) 

CONTINUE 

END  CONDITIONS 

CONTINUE 

IF  (T.LE.0.2)  THEN 
POftM . 744C75-2 . 472E76BT 
ELBE 
P0ft»0. 

IF  <T. LE.0.2. Oft. T.OE. 1.45)  THEN 
P0R-5.91E75 

EL IE IF  <T.GT.0.2.ANB.T.LE.0.5S)  THEN 
POft-5.7 1E75-1 . 23f5E76*< T-0 . 2) 

EL SC IF  <T.OT.O.SS.ANB.T.LT.1.4S>  THEN 
POft-l .374E75-3.268E74«<T-0.55> 

ELSE 

P0ft-1.079E7542.414E74*<T-1.45) 

ENB1F 

PPU,l)»POft 

>1 

I-I 

CALL  FCH(I» J»0»8U»P»Bfft»FFfZETA?CN) 
QP<l,l)»CIHftP(l,l)M<l) 

0PU(l,l)-flP(lfl> 

J-N 

I«NS(N) 

CALL  FCP<I» J»QtQU»P»B»ft»FF»ZETAfCP) 

PP(J,I)-P<J,I) 

QP<JI)<P-PP(J,I)*B(J) 

I-T/TC7 

IP  (I.8E.NTA)  GOTO  79 

TIMT-<I-i>tTC>/TC 

TAU<N,2)«TA(I>B(l-TH)7TH»TAdfI) 

IF  (T.8T.TC)  GOTO  79 
TAU<N,2)«(I-T/TC)«BTX 
GOTO  GO 

TAU<H»2)«TA(NTA) 

CV<H»2)*CVP<H»2)lTAU<Hf 2)S*2/B(H) 

0P<NrNG(H))»“CV(H»2)7«OftT<  <CVCH»2))tB272tCV<H»2)tCP) 
PP<H»NG(H) )*<CP“0P<N»NG<N)) )/B<H) 


C  X»CV<H»2)/2 

C  PP<N*NS(N> )*<-S8RT(X)/B(H)+8QRT(X+B(H)9(CP+QVtSIN(UtT) ) ) 
C  •  /B<H) )t$2 

C  flr<N,HS(M))-CP-PP(H,HS<H))tB(«) 

C  OOMOOtSOPT<PP(HfM8<H)  >/P0) 

QPU<H*N8<H))*QP(H*MS(H) ) 

IF  (HB.NE.H)  THEN 
DO  OS  >H1,HB 
I-NS(J) 

CALL  FCP<I,J,Q,QU,P,B,R,FF,ZETA,CP> 

OP( J»I)*0. 

qpu<:  j»d«o» 

PP<J,I)»CP 

89  CONTINUE 
ENBIF 

C 

C  RESET  VALUES  FOR  NEXT  ITERATION 
C 

DO  90  >1,HB 

DO  90  I«1*NS<J) 

0Q<J,I)-0(J,I) 

00U( JrI)*OU( Jf I) 

0<J,I)-0P<J,I) 

QU< J*I)«QPU(J*I) 

P<J,I)«PP<J,I> 

IF  (P< J*I) «QTiPHAX)  THEN 
PHAX"P< J»I) 

IHAX-I 

JHAX-J 

TNX-T 

ENBIF 

IF  <P<J,I)fPB.LT.PV)  THEN 
PCJ» I)«-PB9FV 
ENOIF 

90  CONTINUE 

IF  <K/IPR*IPR-K.EO.O>  THEN 
GOTO  40 
ELSE 
GOTO  SO 
ENBIF 
99  STOP 

am 

c 

c 

c 

SUBROUTINE  FCH<I,J,Q,QU,P,B,R,FF,ZETA,CN) 

REAL  0<17*190>*P(17*130)tZETA(17)*B(17)*FF(l7) 

REAL  QU< 17*150) *R< 17) 

C 

QS*QU,I)-ZETA<J>S<Q<J,I)H}<J,I+1)) 

PS"P( J*I)-ZETA< J)B(P( J»I)-P( J*I4l) ) 

CN-OS-Rt J)-PS*B<J)-FF<J)iABS(QS>tOS 
RETURN 


•>  a  :  a*t.  I 


SUBROUTINE  FCP( I r J»Q»OU»Pf B» R,FF,ZETAfCP) 

MM.  <Hl7,130>fPU7,130),ZETA<l7>,BU7>,FF<l7> 
ACM.  QU(17fl50)»R<17) 

QA»0( J» I )-ZETA( J)*<0( J»I)-Q( J,I-l ) ) 

PR*P< J»I)-ZETA< J)t<P(J»I)-P( j,l-i) ) 

CP»QR-R ( J ) +PRBB ( J ) -FF ( J ) SABS (QR)tQR 
RETURN 
END 


SUBROUTINE  D£LP< J,A,QPfPP,NS,DP> 

REAL  QP(17,150)»PP<17r150)»A(l7) 

INTEGER  N8<17) 

COHHON/DELTAP/  BO'Bl'B2,B3'QPO,PPO>DPPO 
I-NS(JI) 

DP-DPPO-BOt < ( OP < J , 1 > -QPO > / A < J ) ) t B2-B1 ♦ C PP C J 1 , I > -PRO ) 
•  /(B2+B3*<PP<J1,I>-PP0)> 

RETURN 

END 


SUBROUTINE  VALVE<  J,B,N8,CP,CH,CV,GP,PP,QPU> 

REAL  B<17),CVU7),QPU7,I50>,PP<!7,130>,QPU<17,150> 

INTEGER  N8(17) 

Jl-J-1 

I-NS(Jl) 

X»CP/B<  JD-CH/K  J) 

IF  <X«GE»0)  THEN 

QP<  J1»I)»-CV<  J)*<  1/B<  J)4T/B<  Jl)  HSQRK  <CV<  J)*(  1/B(  J) 
>  f 1/B< Jl) ) )S#2+2tCV( J)fX) 

ELSE 

QP<JI»I)"CV( J)t(l/B( J)tl/B( J1>)-8QRT( (CV( J)t<l/B( J) 

I  M/B<Jl>>>»2-2tCV<J)«X> 

enbzf 

PP<JltX)*<CP-GP< Jlf I))/B<Jl) 

QPU( Jl»I)*OP( Jlrl) 

QP(Jrl)sGP<Jl>I) 

OPU<J,l)-GP<J»l> 

pp(j,n*<Qp(j,n-CH)/B(j) 

RETURN 

END 


SUBROUTINE  BRANCH( J, J2,D,N8,CP,CN,CH1,QP,PP,QPU> 
REAL  B<17),QP<17,150)'QPU<17'1S0)'PP<17'1S0> 
INTEGER  N8<17) 


[ 

“r 

►t 


i-w- 


Cv'*v 


mmm 


non  n  nno  n  non 


I-NS(Jl) 

Pf<Jtl)*<CP-CH-CHl)/<B(J>41<Jl)+I(J2)) 

PP(J2,1>«PP(J,1> 

PP<J1,I>*PP<J,1> 

QP< J» 1 >«PT< J, 1 )SB< J> +CM1 

QP<J2,l)«ff(J2,l)fB<J2)4CM 

QP<J1,I>«CP-PP<J,1)*B<J1> 

QPU( Jrl)«GP< J?l) 

QPU< J2»l)>Bf ( Jf 1) 

0PU<J1,I)«0P<J,1> 

RETURN 

END 


SUBROUTINE  PUHP< J,B,CH,CP,DP,NS,QP,PP,QPU> 
REM.  B< 17) tOf <17»150) »QPU<17,150) »PP( 17, 150) 
INTEGER  NS(17) 


I-NS(Jl) 

PP( J»1>»(CP-CR+B< J1 >t8P)/(B< J)+B< Jl> ) 
PP<J1,I)«PP<J,1)-DP 
QP<J1,I)-CP-B<J1>*PP<J1,I> 
QP<J,1>-QP(J1,I> 

QPU< J,l)aQP( J,l) 

QPU<J1,I)«0P<J,1) 

RETURN 

END 


SUBROUTINE  PULSEt J,B,T,U,NS,GV,CP,CH,Pf ,QP,QPU) 
REAL  B<17)»PP<17,150),QP<17,150),QPU<17,150) 
INTEGER  N8< 17> 

J1«J-1 

I-NS(Jl) 

PP<  J»l)*<CP-CfHJVtABS<SIN(UBT>))/(B<  JHB<  Ji>> 
PP<J1,I)-PP<J,1> 

QP<  J»1)*PP(  J»1)*B(  J)+CN 
QP<JI»I)*CP-PP<J1,I)*B< Jl) 

GPU< Jt l)aGP( Jf 1 ) 

GnHJltX)-GP(Jltl) 


SUBROUTINE  ACCUH< J,B,P,0,CP,CH,DT,NSfPAO,VAO,TOL,PP,QP,OPU> 
REAL  B<17>»T0L<i7)»P<l7»l5O>»0(l7f 150)tPP(l7»l50> »QP< 17, 150) 
REAL  QPU<17,150) »OLDCP<17),OLDCH<17>,PAO<17)»VAO<17) 

INTEGER  NS (17) 


c 


COHHON/OLD/  OLDCH,OLDCP 


Jl-J-1 

Of»ID-<OLDCA<J)fCH>/2 

CPM»-<0LDCP<Jl)+CP>/2 

roi«p(j,n 

pbi-pend 

CALL  PNEN<J'B»P'DT,PAO,VAO,PEND, CANID, CPNID'PN, BIN, BOUT) 
PNt-PN 

PD2-PENM1.  OOOOl 

CALL  PNEN(J,B,P,DT,PA0,VA0,PD2, CANID, CPHID,PN, BIN, QOUT) 

10  A«<PNI-PN)/<PD1-PD2) 

At>PNl-A«PDl 

PEND-Al/U-A) 

PD1-PD2 

PN1-PN 

CALL  PN£U( J,B,P,DT, PAO,VAO, PEND, CANID, CPH1D,PN, BIN, OOUT) 

IF  (ABS(PN-PEND) »OT»TOL( J))  THEN 
PD2*PEND 
GOTO  10 
ENPIF 

PP< J»1)*PEND 

PP< J1»N8( Jl) )*PP( J*l) 

QP(J»1)«Q0UT 

QP<J1,NS<J1))-0IN 

0PU(J,1)-0P<J,1> 

QPU< J1,NS< J1 ) )»QP< J1 » N8< Jl> ) 

OLDCHUXH 
0LDCP<J1>-CP 
RE TUAN 
EM 
C 
C 
C 

SUBROUTINE  PNEH< J,B,P,DT, PAO,VAO, PEND, CANID, CPAID,PN, BIN, BOUT) 
REAL  P(17,1S0),B(17),PA0U7>,UA0(17) 

INTEGER  N8<17) 

C 

Jl-J-1 

PAID»<P<J,l)+PEND>/2 
OIN»CPHID-PAIWB(  Jl) 

BOUT«PHXDt»(J)fCHNID 

VDOT-BIH-flOUT 

PH*  <  J ,  1 HDT  »VDO  T  PPA I  Dt  *  2/  <  P  AO  <  J )  »V  AO  <  J  > ) 
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Appendix  B:  Additional  Boundary  Condition* 


Branch  Connection 

For  the  purpose*  of  this  work,  any  branching  junctions  were 
considered  to  consist  of  one  pip*  flowing  into  two  (Fig.  B-l).  In  this 
type  of  connection  the  continuity  equation  was  used,  a  comon  head  was 
assumed  after  neglecting  losses  as  minor  effects,  and  the  characteristic 
equations  at  the  pipe  ends  were  needed:  Eq.  (2.18)  for  pipe  1,  and  Eq. 
(2.19)  for  pipes  2  and  3  (14:45).  Rewriting  the  characteristic  equations 


Pp-Ppi-Pp2-Pp3 

Qpi-Cp-BPpi 

QP2“cM1+bpP2 

Qp3-C|ffBPp3 


(B.l) 


(B.2) 
( B .  3 ) 
(B.4) 


A  sussaation  of  Eq.  (B.l)  through  (B.4)  provided  a  simple  solution 


for  the  common  pressure,  Pp: 

Pp-(Cp-CHrCM)/2B 

The  flow  in  each  pipe  would  then  be  found  by  us*  of  the  appropriate 


(B.5) 


characteristic  equation.  If  required,  this  method  could  be  expanded  to 
any  other  number  of  pipes  by  the  addition  of  more  characteristic 
equations  to  the  summation  in  Eq.  (B.5). 


Fig.  B-l.  Branch  Connection 


,'t\H  , 


WWW 


Initiatmoui  Pressure  Rise 

The  end  conditions  created  by  e  nonreciproceting  puap  aay  be 
approximated  by  a  luaped  rise  in  pressure  of  aagnitude,  AP ,  as  long  as 
the  diaensions  of  the  puap  are  saall  coapared  to  those  of  the  other 
systea  eleaents  (3:16).  The  end  conditions  take  the  form 


Qpi-Cp-BiPpi  QP2“CM~B2PP2 
QP1"QP2  PP2*PP1+  &P 


(B.5) 


Solving  siaultaneoualy 

PP1"^P“^M“®2  AP)/(Bi+B2)  (B.6) 

Here  again  the  remaining  unknowns  would  be  calculated  directly  from 
the  appropriate  equation.  Note,  however,  that  no  loss  tera  was  used  at 
the  interface.  These  losses  should  be  taken  into  account  when  assigning 
the  functional  fora  of  AP.  One  such  fora  of  AP  was  given  by  Fashbaugh 
and  Streeter  (11:1014)  as 

AP-AP0-B0 ( Qd-Qdo )/A-[B1~(Pu-Pu0)/(B2+(Pu-Puo)B3) ) 2  ( B . 8 ) 

in  which  Qd  is  the  outlet  flow  and  Pu  the  pressure  at  the  inlet  of  the 
puap.  The  constants  AP0,  B0,  and  Qdo  are  constants  determined  from  the 
steady-state  pressure  rise  versus  flow  rate  curve.  While  Puo,  Bj,  B2, 
and  B3  are  constants  obtained  froa  the  steady-state  pressure  rise  versus 
inlet  pressure  curve  (11:1014). 
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